S1: behavioural, conventional analysis with brms
I. S. Plank
2025-09-14
1 Task info and setup
This R Markdown script analyses data from the PAL (probabilistic associative learning) task of the EMBA project. It focuses on the comparison of adults with ADHD (with or without comorbid ASD) and adults without psychiatric diagnoses. The data was preprocessed before being read into this script.
In this task, participants view faces for which they have to judge whether the portrayed emotion is positive or negative. The faces are preceded by a tone which is predictive of the valence of the following face. The visual angle of the faces was 3.46 degrees wide and 4.72 degrees high.
1.1 Plot the development of associations across task
# load the trial order
df.trl = read_csv('../data/PAL_scheme.csv', show_col_types = F) %>%
mutate(
association = ut*100,
prob_HP = prob_HP*100
)
# create rectangles to colour the phases: prevolatile, volatile and postvolatile
rects = data.frame(xstart = c(-Inf, 72, 264),
xend = c(72, 264, Inf),
col = c("a", "b", "c"))
ggplot() +
geom_rect(data = rects, aes(xmin = xstart, xmax = xend,
ymin = -Inf, ymax = Inf,
fill = col), alpha = 0.4) +
scale_fill_manual(values=c(c_mid_highlight, c_green, c_dark_highlight),
guide = "none") +
geom_jitter(data = df.trl,
aes(x = trl, y = association, shape = expected),
colour = "black", alpha = 0.8, size = 1, width = 0, height = 4) +
geom_line(data = df.trl,
aes(x = trl, y = prob_HP),
colour = "black") +
scale_shape_manual(values = c(8, 1, 3)) +
theme_bw() +
scale_y_continuous(breaks=c(16.7, 83.3)) +
labs(x = "trial", y = "percent",
title = "High tone > positive, low tone > negative emotion") +
annotate("text", x = 36, y = 125, label= "prevolatile", size = 4) +
annotate("text", x = 168, y = 125, label= "volatile", size = 4) +
annotate("text", x = 301, y = 125, label= "postvolatile", size = 4) +
theme(legend.position = c(0.9, 0.5), legend.title = element_blank(),
legend.background = element_rect(fill = alpha("white", 0)),
legend.key = element_rect(fill = alpha("white", 0)),
plot.title = element_text(hjust = 0.5), text = element_text(size = 12))
ggsave("plots/PAL_scheme.png",
units = "cm",
width = 27,
height = 9)1.2 Some general settings
# number of simulations
nsim = 250
# set number of iterations and warmup for models
iter = 6000
warm = 1500
# set the seed
set.seed(2468)1.3 Package versions
The following packages are used in this RMarkdown file:
## [1] "R version 4.5.1 (2025-06-13)"## [1] "knitr version 1.50"
## [1] "ggplot2 version 3.5.2"
## [1] "brms version 2.22.0"
## [1] "bridgesampling version 1.1.2"
## [1] "tidyverse version 2.0.0"
## [1] "ggpubr version 0.6.1"
## [1] "ggrain version 0.0.4"
## [1] "bayesplot version 1.13.0"
## [1] "SBC version 0.3.0.9000"
## [1] "rstatix version 0.7.2"
## [1] "officer version 0.6.10"
## [1] "effectsize version 1.0.1"
## [1] "BayesFactor version 0.9.12.4.7"
## [1] "bayestestR version 0.17.0"1.4 Analysis guidelines
We planned to determine the group-level effect subjects following Barr (2013). For each model, experiment specific priors were set based on previous literature or the task (see comments in the code).
We performed prior predictive checks as proposed in Schad, Betancourt and Vasishth (2020) using the SBC package before running this analysis (see [!LINK]). We performed these data-naive evaluations based on the comparison of three groups. To do so, we create 250 simulated datasets where parameters are simulated from the priors. These parameters are used to create one fake dataset. Both the true underlying parameters and the simulated discrimination values are saved. Then, we create graphs showing the prior predictive distribution of the simulated discrimination threshold to check whether our priors fit our general expectations about the data. Next, we perform checks of computational faithfulness and model sensitivity as proposed by Schad, Betancourt and Vasishth (2020) and implemented in the SBC package. We create models for each of the simulated datasets. Last, we calculate performance metrics for each of these models, focusing on the population-level parameters.
1.5 Preparation
First, we load the data and combine it with demographic information including the diagnostic status of the subjects. Then, all predictors are set to sum contrasts.
# check if the data file exists, if yes load it:
if (!file.exists("../data/PAL-ADHD_data.RData")) {
path = "/home/emba/Documents/EMBA"
# get demo info for subjects
df.sub = read_csv(file.path(path, "CentraXX", "EMBA_centraXX.csv"),
show_col_types = F) %>%
mutate(
diagnosis = as.factor(recode(diagnosis, "CTR" = "COMP"))
)
# set the data path
dt.path = file.path(path, "BVET")
dt.explo = file.path(path, "BVET-explo")
# load the preprocessed data: df.tsk
df.tsk = rbind(readRDS(file.path(dt.path, "PAL_tsk.rds")),
readRDS(file.path(dt.explo, "PAL_tsk.rds")))
# load excluded participants (accuracy < 2/3)
exc = c(scan(file.path(dt.path, 'PAL_exc.txt'), what="character", sep=NULL),
scan(file.path(dt.explo, 'PAL_exc.txt'), what="character", sep=NULL))
# merge with group and focus on comparison and autistic adults
df.tsk = merge(df.sub %>% select(subID, diagnosis),
df.tsk, all.y = T) %>%
filter(diagnosis != "ASD" & !(subID %in% exc)) %>%
droplevels() %>%
mutate_if(is.character, as.factor)
# only keep participants included in the study in the subject data frame
df.sub = merge(df.tsk %>% select(subID) %>% distinct(), df.sub, all.x = T) %>%
droplevels()
# anonymise the data based on predetermined file > same for all data
df.recode = read_csv(file.path(path, "BVET", "PID_anonymisation.csv"))
recode = as.character(df.recode$subID)
names(recode) = df.recode$PID
df.tsk$subID = str_replace_all(df.tsk$subID, recode)
# save the data for use with the EMBA_HGF toolbox
df.tsk %>% select(subID, trl, diagnosis, ut, difficulty, acc, rt) %>%
ungroup() %>%
arrange(subID, trl) %>%
write_csv(., file = "../data/PAL-ADHD_data.csv")
# print gender frequencies and compare them across groups
tb.gen = xtabs(~ gender + diagnosis, data = df.sub)
ct.full = contingencyTableBF(tb.gen,
sampleType = "indepMulti",
fixedMargin = "cols")
# get all self-descriptions for DAN group
gen.desc = unique(tolower(df.sub[df.sub$gender == "dan",]$gender_desc))
# create a demographics overview
ls.demo = createDFdemo(df.sub,
c("age", "edu", "BDI_total", "ASRS_total", "RAADS_total", "TAS_total", "iq"),
"diagnosis")
# save it all
save(df.tsk, ls.demo, ct.full, gen.desc, tb.gen,
file = "../data/PAL-ADHD_data.RData")
} else {
load("../data/PAL-ADHD_data.RData")
}
# print the outcome of the two contingency tables for comparison: all participants
ct.full@bayesFactor## bf error time code
## Non-indep. (a=1) -2.156359 0 Fri Sep 5 12:54:03 2025 21ea4a28ddd8# print gender and diagnosis distribution
kable(tb.gen)| ADHD | BOTH | COMP | |
|---|---|---|---|
| dan | 2 | 3 | 0 |
| fem | 9 | 12 | 10 |
| mal | 11 | 7 | 12 |
# drop the neutral condition for the analysis
df.pal = df.tsk %>%
filter(expected != "neutral" & !is.na(rt.cor)) %>% droplevels() %>%
mutate_if(is.character, as.factor) %>%
ungroup()
# calculate variance of reaction times:
# no difficulty due to suboptimal posterior fit
df.var = df.tsk %>%
filter(expected != "neutral") %>% droplevels() %>%
group_by(subID, diagnosis, phase, expected) %>%
summarise(
totaln = n(),
valuen = sum(!is.na(rt.cor)),
rt.var = sd(rt.cor, na.rm = T)
) %>%
mutate(
perc = valuen/totaln
) %>% filter(perc >= 2/3) %>%
mutate_if(is.character, as.factor) %>% ungroup()
# set and print the contrasts
contrasts(df.pal$diagnosis) = contr.sum(3)
contrasts(df.pal$diagnosis)## [,1] [,2]
## ADHD 1 0
## BOTH 0 1
## COMP -1 -1contrasts(df.pal$expected) = contr.sum(2)
contrasts(df.pal$expected)## [,1]
## expected 1
## unexpected -1contrasts(df.pal$phase) = contr.sum(3)
contrasts(df.pal$phase)## [,1] [,2]
## prevolatile 1 0
## volatile 0 1
## postvolatile -1 -1contrasts(df.pal$difficulty) = contr.sum(3)
contrasts(df.pal$difficulty)## [,1] [,2]
## easy 1 0
## medium 0 1
## difficult -1 -1contrasts(df.var$diagnosis) = contr.sum(3)
contrasts(df.var$diagnosis)## [,1] [,2]
## ADHD 1 0
## BOTH 0 1
## COMP -1 -1contrasts(df.var$expected) = contr.sum(2)
contrasts(df.var$expected)## [,1]
## expected 1
## unexpected -1contrasts(df.var$phase) = contr.sum(3)
contrasts(df.var$phase)## [,1] [,2]
## prevolatile 1 0
## volatile 0 1
## postvolatile -1 -1# print final group comparisons for the paper into a Word document
read_docx() %>%
body_add_table(ls.demo$df.demo %>% arrange(measurement) %>%
mutate(bf.log =
if_else(
bf.log > log(3),
sprintf("%.3f*", bf.log),
sprintf("%.3f", bf.log)))) %>%
print(target = "PAL-ADHD_demo.docx")
read_docx() %>%
body_add_table(ls.demo$df.post %>% arrange(measurement) %>%
mutate(bf.log =
if_else(
bf.log > log(3),
sprintf("%.3f*", bf.log),
sprintf("%.3f", bf.log)))) %>%
print(target = "PAL-ADHD_post.docx")2 Reaction time variance
In the preregistration, we noted the following population-level effects for the model of the reaction time variances: group, expectancy, phase and difficulty. However, the posterior fit for this model was suboptimal, therefore, we dropped the predictor difficulty.
2.1 Model setup
# figure out slopes for subjects
kable(head(df.var %>% count(subID, expected)))| subID | expected | n |
|---|---|---|
| 1 | expected | 3 |
| 1 | unexpected | 3 |
| 10 | expected | 3 |
| 10 | unexpected | 3 |
| 11 | expected | 3 |
| 11 | unexpected | 3 |
kable(head(df.var %>% count(subID, phase)))| subID | phase | n |
|---|---|---|
| 1 | prevolatile | 2 |
| 1 | volatile | 2 |
| 1 | postvolatile | 2 |
| 10 | prevolatile | 2 |
| 10 | volatile | 2 |
| 10 | postvolatile | 2 |
kable(head(df.var %>% count(subID, expected, phase)))| subID | expected | phase | n |
|---|---|---|---|
| 1 | expected | prevolatile | 1 |
| 1 | expected | volatile | 1 |
| 1 | expected | postvolatile | 1 |
| 1 | unexpected | prevolatile | 1 |
| 1 | unexpected | volatile | 1 |
| 1 | unexpected | postvolatile | 1 |
# model formula
f.var = brms::bf(rt.var ~ diagnosis * expected * phase +
(expected + phase | subID)
)
# set weakly informative priors
priors = c(
prior(normal(4.5, 0.50), class = Intercept),
prior(normal(0, 0.50), class = sigma),
prior(normal(0.15, 0.15), class = sd),
prior(normal(0, 0.25), class = b),
prior(lkj(2), class = cor)
)2.2 Posterior predictive checks
As the next step, we fit the model, check whether there are divergence or rhat issues, and then check whether the chains have converged.
# fit the final model
m.var = brm(f.var, seed = 4646,
df.var, prior = priors,
family = lognormal,
iter = iter, warmup = warm,
backend = "cmdstanr", threads = threading(8),
file = file.path(brms_dir, "m_var"),
save_pars = save_pars(all = TRUE)
)
rstan::check_hmc_diagnostics(m.var$fit)##
## Divergences:## 0 of 18000 iterations ended with a divergence.##
## Tree depth:## 0 of 18000 iterations saturated the maximum tree depth of 10.##
## Energy:## E-BFMI indicated no pathological behavior.# check that rhats are below 1.01
sum(brms::rhat(m.var) >= 1.01, na.rm = T)## [1] 0# check the trace plots
post.draws = as_draws_df(m.var)
mcmc_trace(post.draws, regex_pars = "^b_",
facet_args = list(ncol = 4)) +
scale_x_continuous(breaks=scales::pretty_breaks(n = 3)) +
scale_y_continuous(breaks=scales::pretty_breaks(n = 3))## Scale for x is already present.
## Adding another scale for x, which will replace the existing scale.
This model has no pathological behaviour with E-BFMI, no divergent sample and no rhats that are higher or equal to 1.01. Therefore, we go ahead and perform our posterior predictive checks.
# get posterior predictions
post.pred = posterior_predict(m.var, ndraws = nsim)
# check the fit of the predicted data compared to the real data
p1 = pp_check(m.var, ndraws = nsim) +
theme_bw() + theme(legend.position = "none") + xlim(0, 1000)
# distributions of means compared to the real values per group
p2 = ppc_stat_grouped(df.var$rt.var, post.pred, df.var$diagnosis) +
theme_bw() + theme(legend.position = "none")
# distributions of means compared to the real values per phase
p3 = ppc_stat_grouped(df.var$rt.var, post.pred, df.var$phase) +
theme_bw() + theme(legend.position = "none")
# distributions of means compared to the real values per expected
p4 = ppc_stat_grouped(df.var$rt.var, post.pred, df.var$expected) +
theme_bw() + theme(legend.position = "none")
p = ggarrange(p1, p2, p3, p4,
nrow = 4, ncol = 1, labels = "AUTO")
annotate_figure(p, top = text_grob("Posterior predictive checks: RT variance",
face = "bold", size = 14))
The simulated data based on the model (light blue) fits well with the real data (dark blue), both for the overall distribution (A) and the aggregated scores for the diagnostic groups (B).
2.3 Inferences
Now that we are convinced that we can trust our model, we have a look at its estimate and use the hypothesis function to assess our hypotheses and perform explorative tests.
# print a summary
summary(m.var)## Family: lognormal
## Links: mu = identity; sigma = identity
## Formula: rt.var ~ diagnosis * expected * phase + (expected + phase | subID)
## Data: df.var (Number of observations: 387)
## Draws: 4 chains, each with iter = 6000; warmup = 1500; thin = 1;
## total post-warmup draws = 18000
##
## Multilevel Hyperparameters:
## ~subID (Number of levels: 66)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## sd(Intercept) 0.36 0.03 0.30 0.43 1.00 3864
## sd(expected1) 0.02 0.01 0.00 0.05 1.00 4508
## sd(phase1) 0.09 0.02 0.05 0.12 1.00 6414
## sd(phase2) 0.03 0.02 0.00 0.08 1.00 4058
## cor(Intercept,expected1) -0.20 0.31 -0.74 0.47 1.00 22631
## cor(Intercept,phase1) -0.23 0.18 -0.56 0.13 1.00 17311
## cor(expected1,phase1) 0.26 0.33 -0.48 0.81 1.00 2930
## cor(Intercept,phase2) 0.04 0.30 -0.56 0.62 1.00 22936
## cor(expected1,phase2) 0.19 0.37 -0.59 0.81 1.00 6950
## cor(phase1,phase2) 0.08 0.34 -0.58 0.70 1.00 17998
## Tail_ESS
## sd(Intercept) 7484
## sd(expected1) 5809
## sd(phase1) 6234
## sd(phase2) 6446
## cor(Intercept,expected1) 12093
## cor(Intercept,phase1) 13589
## cor(expected1,phase1) 5561
## cor(Intercept,phase2) 12392
## cor(expected1,phase2) 11325
## cor(phase1,phase2) 14323
##
## Regression Coefficients:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## Intercept 4.71 0.04 4.62 4.79 1.00 2209
## diagnosis1 0.09 0.06 -0.03 0.21 1.00 2291
## diagnosis2 0.03 0.06 -0.09 0.15 1.00 2444
## expected1 0.06 0.01 0.04 0.08 1.00 30294
## phase1 0.02 0.02 -0.01 0.05 1.00 15167
## phase2 0.03 0.01 0.00 0.06 1.00 26066
## diagnosis1:expected1 -0.01 0.01 -0.04 0.02 1.00 22496
## diagnosis2:expected1 -0.00 0.01 -0.03 0.02 1.00 22150
## diagnosis1:phase1 0.00 0.02 -0.05 0.05 1.00 13212
## diagnosis2:phase1 0.05 0.02 -0.00 0.10 1.00 13450
## diagnosis1:phase2 -0.03 0.02 -0.07 0.01 1.00 19605
## diagnosis2:phase2 -0.02 0.02 -0.06 0.02 1.00 20126
## expected1:phase1 0.02 0.01 -0.00 0.05 1.00 27027
## expected1:phase2 -0.02 0.01 -0.05 0.00 1.00 25645
## diagnosis1:expected1:phase1 -0.01 0.02 -0.05 0.03 1.00 19059
## diagnosis2:expected1:phase1 -0.01 0.02 -0.05 0.02 1.00 19306
## diagnosis1:expected1:phase2 0.04 0.02 -0.00 0.07 1.00 19797
## diagnosis2:expected1:phase2 0.00 0.02 -0.03 0.04 1.00 18557
## Tail_ESS
## Intercept 5022
## diagnosis1 4509
## diagnosis2 4865
## expected1 13918
## phase1 13959
## phase2 13147
## diagnosis1:expected1 14370
## diagnosis2:expected1 14145
## diagnosis1:phase1 13651
## diagnosis2:phase1 13448
## diagnosis1:phase2 13824
## diagnosis2:phase2 14268
## expected1:phase1 15127
## expected1:phase2 15707
## diagnosis1:expected1:phase1 15059
## diagnosis2:expected1:phase1 13226
## diagnosis1:expected1:phase2 14941
## diagnosis2:expected1:phase2 14865
##
## Further Distributional Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 0.18 0.01 0.17 0.20 1.00 6545 10033
##
## Draws were sampled using sample(hmc). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).# get the estimates and compute group comparisons
df.m.var = post.draws %>%
select(starts_with("b_")) %>%
mutate(
b_COMP = - b_diagnosis1 - b_diagnosis2,
b_postvolatile = - b_phase1 - b_phase2,
BOTH = b_Intercept + b_diagnosis2,
ADHD = b_Intercept + b_diagnosis1,
COMP = b_Intercept + b_COMP,
`ADHD-BOTH` = ADHD - BOTH,
`ADHD-COMP` = ADHD - COMP,
`BOTH-COMP` = BOTH - COMP
)
# plot the posterior distributions
df.m.var %>%
select(starts_with("b_")) %>%
pivot_longer(cols = starts_with("b_"), names_to = "coef", values_to = "estimate") %>%
subset(!startsWith(coef, "b_Int")) %>%
mutate(
coef = substr(coef, 3, nchar(coef)),
coef = str_replace_all(coef, ":", " x "),
coef = str_replace_all(coef, "diagnosis1", "ADHD"),
coef = str_replace_all(coef, "diagnosis2", "ADHD+ASD"),
coef = str_replace_all(coef, "expected1", "expected"),
coef = str_replace_all(coef, "expected2", "unexpected"),
coef = str_replace_all(coef, "phase1", "prevolatile"),
coef = str_replace_all(coef, "phase2", "volatile"),
coef = fct_reorder(coef, desc(estimate))
) %>%
group_by(coef) %>%
mutate(
cred = case_when(
(mean(estimate) < 0 & quantile(estimate, probs = 0.975) < 0) |
(mean(estimate) > 0 & quantile(estimate, probs = 0.025) > 0) ~ "credible",
T ~ "not credible"
)
) %>% ungroup() %>%
ggplot(aes(x = estimate, y = coef, fill = cred)) +
geom_vline(xintercept = 0, linetype = 'dashed') +
ggdist::stat_halfeye(alpha = 0.7) + ylab(NULL) + theme_bw() +
scale_fill_manual(values = c(c_dark, c_light)) + theme(legend.position = "none")
# H1a: COMP < ADHD
h1a = hypothesis(m.var, "0 < 2*diagnosis1 + diagnosis2")
h1a## Hypothesis Tests for class b:
## Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio
## 1 (0)-(2*diagnosis1... < 0 -0.21 0.11 -0.39 -0.03 37.54
## Post.Prob Star
## 1 0.97 *
## ---
## 'CI': 90%-CI for one-sided and 95%-CI for two-sided hypotheses.
## '*': For one-sided hypotheses, the posterior probability exceeds 95%;
## for two-sided hypotheses, the value tested against lies outside the 95%-CI.
## Posterior probabilities of point hypotheses assume equal prior probabilities.## Exploration
# COMP < BOTH
e1 = hypothesis(m.var, "0 < diagnosis1 + 2*diagnosis2", alpha = 0.025)
e1## Hypothesis Tests for class b:
## Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio
## 1 (0)-(diagnosis1+2... < 0 -0.15 0.11 -0.37 0.06 11.82
## Post.Prob Star
## 1 0.92
## ---
## 'CI': 95%-CI for one-sided and 97.5%-CI for two-sided hypotheses.
## '*': For one-sided hypotheses, the posterior probability exceeds 97.5%;
## for two-sided hypotheses, the value tested against lies outside the 97.5%-CI.
## Posterior probabilities of point hypotheses assume equal prior probabilities.# ADHD vs. BOTH
e2 = hypothesis(m.var, "diagnosis1 > diagnosis2", alpha = 0.025)
e2## Hypothesis Tests for class b:
## Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio
## 1 (diagnosis1)-(dia... > 0 0.06 0.11 -0.15 0.27 2.53
## Post.Prob Star
## 1 0.72
## ---
## 'CI': 95%-CI for one-sided and 97.5%-CI for two-sided hypotheses.
## '*': For one-sided hypotheses, the posterior probability exceeds 97.5%;
## for two-sided hypotheses, the value tested against lies outside the 97.5%-CI.
## Posterior probabilities of point hypotheses assume equal prior probabilities.## Exploration of possible transition effects
## create design matrix to figure out how to set comparisons for hypotheses
df.des = cbind(df.var,
model.matrix(~ diagnosis * expected * phase,
data = df.var)) %>%
select(-subID, -rt.var, -perc, -totaln, -valuen) %>% distinct()
# COMP(vol-pre) > ADHD(vol-pre)
as.data.frame(t(df.des %>%
ungroup() %>%
filter((diagnosis == "ADHD" | diagnosis == "COMP") &
((phase == "volatile") | (phase == "prevolatile"))) %>%
group_by(diagnosis, phase) %>%
summarise(across(where(is.numeric), ~ mean(.x))) %>%
group_by(diagnosis) %>%
summarise(across(where(is.numeric), ~ diff(.x))) %>%
select(where(is.numeric)) %>%
map_df(~ diff(.x)))) %>%
filter(V1 != 0)## V1
## diagnosis1:phase1 2
## diagnosis2:phase1 1
## diagnosis1:phase2 -2
## diagnosis2:phase2 -1e3 = hypothesis(m.var, "0 < 2*diagnosis1:phase1 + diagnosis2:phase1 +
- (2*diagnosis1:phase2 + diagnosis2:phase2)",
alpha = 0.025)
e3## Hypothesis Tests for class b:
## Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio
## 1 (0)-(2*diagnosis1... < 0 -0.13 0.06 -0.26 -0.01 46.62
## Post.Prob Star
## 1 0.98 *
## ---
## 'CI': 95%-CI for one-sided and 97.5%-CI for two-sided hypotheses.
## '*': For one-sided hypotheses, the posterior probability exceeds 97.5%;
## for two-sided hypotheses, the value tested against lies outside the 97.5%-CI.
## Posterior probabilities of point hypotheses assume equal prior probabilities.# COMP(vol-pre) > BOTH(vol-pre)
e4 = hypothesis(m.var, "0 < diagnosis1:phase1 + 2*diagnosis2:phase1 +
- (diagnosis1:phase2 + 2*diagnosis2:phase2)",
alpha = 0.025)
e4## Hypothesis Tests for class b:
## Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio
## 1 (0)-(diagnosis1:p... < 0 -0.17 0.06 -0.29 -0.04 201.25
## Post.Prob Star
## 1 1 *
## ---
## 'CI': 95%-CI for one-sided and 97.5%-CI for two-sided hypotheses.
## '*': For one-sided hypotheses, the posterior probability exceeds 97.5%;
## for two-sided hypotheses, the value tested against lies outside the 97.5%-CI.
## Posterior probabilities of point hypotheses assume equal prior probabilities.# COMP(post-vol) < ADHD(post-vol)
as.data.frame(t(df.des %>%
ungroup() %>%
filter((diagnosis == "ADHD" | diagnosis == "COMP") &
((phase == "volatile") | (phase == "postvolatile"))) %>%
group_by(diagnosis, phase) %>%
summarise(across(where(is.numeric), ~ mean(.x))) %>%
group_by(diagnosis) %>%
summarise(across(where(is.numeric), ~ diff(.x))) %>% # post - volatile
select(where(is.numeric)) %>%
map_df(~ diff(.x)))) %>% # COMP - ADHD
filter(V1 != 0) ## V1
## diagnosis1:phase1 2
## diagnosis2:phase1 1
## diagnosis1:phase2 4
## diagnosis2:phase2 2e5 = hypothesis(m.var, "0 > 2*diagnosis1:phase1 + diagnosis2:phase1 +
+ 2*(2*diagnosis1:phase2 + diagnosis2:phase2)",
alpha = 0.025)
e5## Hypothesis Tests for class b:
## Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio
## 1 (0)-(2*diagnosis1... > 0 0.12 0.07 -0.02 0.25 22.9
## Post.Prob Star
## 1 0.96
## ---
## 'CI': 95%-CI for one-sided and 97.5%-CI for two-sided hypotheses.
## '*': For one-sided hypotheses, the posterior probability exceeds 97.5%;
## for two-sided hypotheses, the value tested against lies outside the 97.5%-CI.
## Posterior probabilities of point hypotheses assume equal prior probabilities.# COMP(post-vol) < BOTH(post-vol)
e6 = hypothesis(m.var, "0 > diagnosis1:phase1 + 2*(diagnosis2:phase1 +
+ diagnosis1:phase2 + 2*diagnosis2:phase2)",
alpha = 0.025)
e6## Hypothesis Tests for class b:
## Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio
## 1 (0)-(diagnosis1:p... > 0 0.05 0.07 -0.08 0.18 3.18
## Post.Prob Star
## 1 0.76
## ---
## 'CI': 95%-CI for one-sided and 97.5%-CI for two-sided hypotheses.
## '*': For one-sided hypotheses, the posterior probability exceeds 97.5%;
## for two-sided hypotheses, the value tested against lies outside the 97.5%-CI.
## Posterior probabilities of point hypotheses assume equal prior probabilities.## extract predicted differences in ms instead of log data
df.new = df.var %>%
select(diagnosis, phase, expected) %>%
distinct() %>%
mutate(
condition = paste(diagnosis, phase, expected, sep = "_")
)
df.ms = as.data.frame(
fitted(m.var, summary = F,
newdata = df.new %>% select(diagnosis, phase, expected),
re_formula = NA))
colnames(df.ms) = df.new$condition
# calculate our difference columns
df.ms = df.ms %>%
mutate(
COMP = rowMeans(across(matches("COMP_.*"))),
ADHD = rowMeans(across(matches("ADHD_.*"))),
BOTH = rowMeans(across(matches("BOTH_.*"))),
h1a_COMPvADHD = ADHD - COMP,
e_COMPvBOTH = BOTH - COMP,
e_BOTHvADHD = BOTH - ADHD,
COMPvol = rowMeans(across(matches("COMP_volatile.*"))),
COMPpre = rowMeans(across(matches("COMP_prevolatile.*"))),
COMPpost = rowMeans(across(matches("COMP_postvolatile.*"))),
ADHDvol = rowMeans(across(matches("ADHD_volatile.*"))),
ADHDpre = rowMeans(across(matches("ADHD_prevolatile.*"))),
ADHDpost = rowMeans(across(matches("ADHD_postvolatile.*"))),
BOTHvol = rowMeans(across(matches("BOTH_volatile.*"))),
BOTHpre = rowMeans(across(matches("BOTH_prevolatile.*"))),
BOTHpost = rowMeans(across(matches("BOTH_postvolatile.*"))),
e_volvpreCOMPvADHD = (COMPvol - COMPpre) - (ADHDvol - ADHDpre),
e_postvvolCOMPvADHD = (COMPpost - COMPvol) - (ADHDpost - ADHDvol),
e_volvpreCOMPvBOTH = (COMPvol - COMPpre) - (BOTHvol - BOTHpre),
e_postvvolCOMPvBOTH = (COMPpost - COMPvol) - (BOTHpost - BOTHvol),
e_volvpreBOTHvADHD = (BOTHvol - BOTHpre) - (ADHDvol - ADHDpre),
e_postvvolBOTHvADHD = (BOTHpost - BOTHvol) - (ADHDpost - ADHDvol)
)
kable(df.ms %>%
select(starts_with("e") | starts_with("h")) %>%
summarise_all(list(lower.ci = lower_ci, mean = mean, upper.ci = upper_ci)) %>%
t %>%
as.data.frame %>%
rownames_to_column() %>%
separate(rowname, into = c("id", "comparison", "stat"), sep = "_") %>%
pivot_wider(names_from = stat, values_from = V1) %>%
arrange(comparison), digits = 2)| id | comparison | lower.ci | mean | upper.ci |
|---|---|---|---|---|
| e | BOTHvADHD | -32.30 | -7.10 | 18.36 |
| h1a | COMPvADHD | -0.58 | 23.43 | 47.23 |
| e | COMPvBOTH | -6.40 | 16.33 | 40.09 |
| e | postvvolBOTHvADHD | -22.91 | -7.22 | 8.29 |
| e | postvvolCOMPvADHD | -26.07 | -10.90 | 3.92 |
| e | postvvolCOMPvBOTH | -17.74 | -3.68 | 10.42 |
| e | volvpreBOTHvADHD | -20.00 | -4.50 | 10.95 |
| e | volvpreCOMPvADHD | -1.24 | 13.23 | 27.92 |
| e | volvpreCOMPvBOTH | 3.57 | 17.72 | 31.92 |
# calculate effect sizes
df.effect = post.draws %>%
mutate(across(starts_with("sd")|starts_with("sigma"), ~.^2)) %>%
mutate(
sumvar = sqrt(rowSums(select(., starts_with("sd")|starts_with("sigma")))),
h1a = (2*b_diagnosis1 + b_diagnosis2) / sumvar,
e1 = (b_diagnosis1 + 2*b_diagnosis2) / sumvar,
e2 = (b_diagnosis1 - b_diagnosis2) / sumvar,
e3 = (`b_diagnosis1:phase1` + 2*`b_diagnosis2:phase1` -
`b_diagnosis1:phase2` - 2*`b_diagnosis2:phase2`) / sumvar,
e4 = (2*`b_diagnosis1:phase1` + `b_diagnosis2:phase1` -
2*`b_diagnosis1:phase2` - `b_diagnosis2:phase2`) / sumvar
)
kable(df.effect %>% select(starts_with("e")|starts_with("h")) %>%
pivot_longer(cols = everything(), values_to = "estimate") %>%
group_by(name) %>%
summarise(
ci.lo = lower_ci(estimate),
mean = mean(estimate),
ci.hi = upper_ci(estimate),
interpret = interpret_cohens_d(mean)
), digits = 3
)| name | ci.lo | mean | ci.hi | interpret |
|---|---|---|---|---|
| e1 | -0.137 | 0.366 | 0.882 | small |
| e2 | -0.356 | 0.147 | 0.642 | very small |
| e3 | 0.099 | 0.401 | 0.711 | small |
| e4 | 0.012 | 0.313 | 0.624 | small |
| h1a | -0.006 | 0.513 | 1.013 | medium |
h1a: estimate = -0.21 [-0.39, -0.03], posterior probability = 97.41%
e1 COMP < BOTH: estimate = -0.15 [-0.37, 0.06], posterior probability = 92.2%
e3 COMP(vol-pre) > ADHD(vol-pre): estimate = -0.13 [-0.26, -0.01], posterior probability = 97.9%
e4 COMP(vol-pre) > BOTH(vol-pre): estimate = -0.17 [-0.29, -0.04], posterior probability = 99.51%
2.4 Plots
# rain cloud plot
df.var %>%
mutate(
group = case_when(
diagnosis == "BOTH" ~ "ADHD+ASD",
T ~ diagnosis
)
) %>%
ggplot(aes(expected, rt.var, fill = group, colour = group)) + #
geom_rain(rain.side = 'r',
boxplot.args = list(color = "black", outlier.shape = NA, show.legend = FALSE, alpha = .8),
violin.args = list(color = "black", outlier.shape = NA, alpha = .8),
boxplot.args.pos = list(
position = ggpp::position_dodgenudge(x = 0, width = 0.3), width = 0.3
),
point.args = list(show_guide = FALSE, alpha = .5),
violin.args.pos = list(
width = 0.6, position = position_nudge(x = 0.16)),
point.args.pos = list(position = ggpp::position_dodgenudge(x = -0.25, width = 0.1))) +
scale_fill_manual(values = col.grp) +
scale_color_manual(values = col.grp) +
#ylim(0, 1) +
labs(title = "", x = "", y = "variance (ms)") +
facet_wrap(. ~ phase) +
theme_bw() +
theme(legend.position = "bottom", plot.title = element_blank(),
legend.direction = "horizontal", text = element_text(size = 15),
legend.title = element_blank())
ggsave("plots/FigRTvar.svg", units = "cm", width = 27, height = 9)
# plot transition effects specifically
df.var %>%
group_by(diagnosis, phase) %>%
summarise(
var.mn = mean(rt.var),
var.se = sd(rt.var) / sqrt(n())
) %>%
ggplot(aes(y = var.mn, x = phase,
group = diagnosis, colour = diagnosis)) +
geom_line(position = position_dodge(0.4), linewidth = 1) +
geom_errorbar(aes(ymin = var.mn - var.se,
ymax = var.mn + var.se), linewidth = 1, width = 0.5,
position = position_dodge(0.4)) +
geom_point(position = position_dodge(0.4), size = 5) +
labs (x = "phase", y = "rt var (ms)") +
scale_fill_manual(values = col.grp) +
scale_color_manual(values = col.grp) +
theme_bw() +
theme(legend.position = "bottom", plot.title = element_text(hjust = 0.5))
3 Reaction times in correct trials
In the preregistration, we noted the following population-level effects for the model of the reaction time variances: group, expectancy, phase and difficulty; as well as the group-level predictores subject and trials.
# figure out slopes for subjects
kable(head(df.pal %>% count(subID, expected)))| subID | expected | n |
|---|---|---|
| 1 | expected | 224 |
| 1 | unexpected | 45 |
| 10 | expected | 225 |
| 10 | unexpected | 43 |
| 11 | expected | 225 |
| 11 | unexpected | 46 |
kable(head(df.pal %>% count(subID, phase)))| subID | phase | n |
|---|---|---|
| 1 | prevolatile | 66 |
| 1 | volatile | 133 |
| 1 | postvolatile | 70 |
| 10 | prevolatile | 67 |
| 10 | volatile | 139 |
| 10 | postvolatile | 62 |
kable(head(df.pal %>% count(subID, difficulty)))| subID | difficulty | n |
|---|---|---|
| 1 | easy | 95 |
| 1 | medium | 94 |
| 1 | difficult | 80 |
| 10 | easy | 91 |
| 10 | medium | 89 |
| 10 | difficult | 88 |
kable(head(df.pal %>% count(subID, expected, phase)))| subID | expected | phase | n |
|---|---|---|---|
| 1 | expected | prevolatile | 55 |
| 1 | expected | volatile | 111 |
| 1 | expected | postvolatile | 58 |
| 1 | unexpected | prevolatile | 11 |
| 1 | unexpected | volatile | 22 |
| 1 | unexpected | postvolatile | 12 |
kable(head(df.pal %>% count(subID, expected, difficulty)))| subID | expected | difficulty | n |
|---|---|---|---|
| 1 | expected | easy | 79 |
| 1 | expected | medium | 78 |
| 1 | expected | difficult | 67 |
| 1 | unexpected | easy | 16 |
| 1 | unexpected | medium | 16 |
| 1 | unexpected | difficult | 13 |
kable(head(df.pal %>% count(subID, phase, difficulty)))| subID | phase | difficulty | n |
|---|---|---|---|
| 1 | prevolatile | easy | 23 |
| 1 | prevolatile | medium | 23 |
| 1 | prevolatile | difficult | 20 |
| 1 | volatile | easy | 48 |
| 1 | volatile | medium | 47 |
| 1 | volatile | difficult | 38 |
kable(head(df.pal %>% count(subID, expected, phase, difficulty)))| subID | expected | phase | difficulty | n |
|---|---|---|---|---|
| 1 | expected | prevolatile | easy | 19 |
| 1 | expected | prevolatile | medium | 19 |
| 1 | expected | prevolatile | difficult | 17 |
| 1 | expected | volatile | easy | 40 |
| 1 | expected | volatile | medium | 39 |
| 1 | expected | volatile | difficult | 32 |
# figure out slopes for trls
kable(head(df.pal %>% count(trl, diagnosis)))| trl | diagnosis | n |
|---|---|---|
| 1 | ADHD | 14 |
| 1 | BOTH | 20 |
| 1 | COMP | 16 |
| 2 | ADHD | 18 |
| 2 | BOTH | 18 |
| 2 | COMP | 21 |
# set the formula
f.pal = brms::bf(rt.cor ~ diagnosis * expected * phase * difficulty +
(expected * phase * difficulty | subID) + (diagnosis | trl))
# set informed priors based on previous results
priors = c(
# informative priors based Lawson et al. and Schad, Betancourt & Vasishth (2019)
prior(normal(6.0, 0.3), class = Intercept),
prior(normal(0.0, 0.5), class = sigma),
prior(normal(0, 0.1), class = sd),
prior(lkj(2), class = cor),
prior(normal(100, 100.0), class = ndt),
prior(normal(0.00, 0.04), class = b)
)3.1 Posterior predictive checks
As the next step, we fit the model, check whether there are divergence or rhat issues, and then check whether the chains have converged.
# fit the final model
m.pal = brm(f.pal, seed = 4466,
df.pal, prior = priors,
family = shifted_lognormal,
iter = iter, warmup = warm,
backend = "cmdstanr", threads = threading(8),
file = file.path(brms_dir, "m_pal"),
save_pars = save_pars(all = TRUE)
)
rstan::check_hmc_diagnostics(m.pal$fit)##
## Divergences:## 0 of 18000 iterations ended with a divergence.##
## Tree depth:## 0 of 18000 iterations saturated the maximum tree depth of 10.##
## Energy:## E-BFMI indicated no pathological behavior.# check that rhats are below 1.01
sum(brms::rhat(m.pal) >= 1.01, na.rm = T)## [1] 0# check the trace plots
post.draws = as_draws_df(m.pal)
mcmc_trace(post.draws, regex_pars = "^b_",
facet_args = list(ncol = 6)) +
scale_x_continuous(breaks=scales::pretty_breaks(n = 3)) +
scale_y_continuous(breaks=scales::pretty_breaks(n = 3))## Scale for x is already present.
## Adding another scale for x, which will replace the existing scale.
This model has no pathological behaviour with E-BFMI, no divergent sample and no rhats that are higher or equal to 1.01. Therefore, we go ahead and perform our posterior predictive checks.
# get posterior predictions
post.pred = posterior_predict(m.pal, ndraws = nsim)
# check the fit of the predicted data compared to the real data
p1 = pp_check(m.pal, ndraws = nsim) +
theme_bw() + theme(legend.position = "none") + xlim(0, 1500)
# distributions of means compared to the real values per group or conditions
p2 = ppc_stat_grouped(df.pal$rt.cor,
post.pred,
df.pal$diagnosis) +
theme_bw() + theme(legend.position = "none")
p3 = ppc_stat_grouped(df.pal$rt.cor,
post.pred,
df.pal$expected) +
theme_bw() + theme(legend.position = "none")
p4 = ppc_stat_grouped(df.pal$rt.cor,
post.pred,
df.pal$phase) +
theme_bw() + theme(legend.position = "none")
p5 = ppc_stat_grouped(df.pal$rt.cor,
post.pred,
df.pal$difficulty) +
theme_bw() + theme(legend.position = "none")
p = ggarrange(p1, p2, p3, p4, p5,
nrow = 5, ncol = 1, labels = "AUTO")
annotate_figure(p, top = text_grob("Posterior predictive checks: reaction times",
face = "bold", size = 14))
The simulated data based on the model fits well with the real data, although there are slight deviations specifically for the predictor expectedness. However, we judge this to be acceptable.
3.2 Inferences
Now that we are convinced that we can trust our model, we have a look at its estimate and use the hypothesis function to assess our hypotheses and perform explorative tests.
# print a summary
summary(m.pal)## Family: shifted_lognormal
## Links: mu = identity; sigma = identity; ndt = identity
## Formula: rt.cor ~ diagnosis * expected * phase * difficulty + (expected * phase * difficulty | subID) + (diagnosis | trl)
## Data: df.pal (Number of observations: 16886)
## Draws: 4 chains, each with iter = 6000; warmup = 1500; thin = 1;
## total post-warmup draws = 18000
##
## Multilevel Hyperparameters:
## ~subID (Number of levels: 66)
## Estimate
## sd(Intercept) 0.19
## sd(expected1) 0.02
## sd(phase1) 0.06
## sd(phase2) 0.02
## sd(difficulty1) 0.01
## sd(difficulty2) 0.01
## sd(expected1:phase1) 0.02
## sd(expected1:phase2) 0.01
## sd(expected1:difficulty1) 0.00
## sd(expected1:difficulty2) 0.01
## sd(phase1:difficulty1) 0.01
## sd(phase2:difficulty1) 0.01
## sd(phase1:difficulty2) 0.01
## sd(phase2:difficulty2) 0.01
## sd(expected1:phase1:difficulty1) 0.01
## sd(expected1:phase2:difficulty1) 0.01
## sd(expected1:phase1:difficulty2) 0.01
## sd(expected1:phase2:difficulty2) 0.01
## cor(Intercept,expected1) 0.23
## cor(Intercept,phase1) 0.08
## cor(expected1,phase1) 0.22
## cor(Intercept,phase2) 0.09
## cor(expected1,phase2) 0.04
## cor(phase1,phase2) -0.25
## cor(Intercept,difficulty1) -0.25
## cor(expected1,difficulty1) -0.08
## cor(phase1,difficulty1) -0.17
## cor(phase2,difficulty1) -0.07
## cor(Intercept,difficulty2) -0.04
## cor(expected1,difficulty2) -0.11
## cor(phase1,difficulty2) -0.17
## cor(phase2,difficulty2) 0.09
## cor(difficulty1,difficulty2) 0.02
## cor(Intercept,expected1:phase1) 0.03
## cor(expected1,expected1:phase1) 0.17
## cor(phase1,expected1:phase1) 0.09
## cor(phase2,expected1:phase1) -0.17
## cor(difficulty1,expected1:phase1) 0.06
## cor(difficulty2,expected1:phase1) -0.02
## cor(Intercept,expected1:phase2) 0.21
## cor(expected1,expected1:phase2) 0.07
## cor(phase1,expected1:phase2) -0.05
## cor(phase2,expected1:phase2) 0.04
## cor(difficulty1,expected1:phase2) 0.05
## cor(difficulty2,expected1:phase2) 0.13
## cor(expected1:phase1,expected1:phase2) 0.01
## cor(Intercept,expected1:difficulty1) -0.07
## cor(expected1,expected1:difficulty1) 0.03
## cor(phase1,expected1:difficulty1) 0.03
## cor(phase2,expected1:difficulty1) 0.00
## cor(difficulty1,expected1:difficulty1) -0.01
## cor(difficulty2,expected1:difficulty1) -0.00
## cor(expected1:phase1,expected1:difficulty1) 0.04
## cor(expected1:phase2,expected1:difficulty1) 0.01
## cor(Intercept,expected1:difficulty2) 0.01
## cor(expected1,expected1:difficulty2) 0.05
## cor(phase1,expected1:difficulty2) 0.10
## cor(phase2,expected1:difficulty2) 0.07
## cor(difficulty1,expected1:difficulty2) 0.02
## cor(difficulty2,expected1:difficulty2) -0.01
## cor(expected1:phase1,expected1:difficulty2) 0.06
## cor(expected1:phase2,expected1:difficulty2) 0.07
## cor(expected1:difficulty1,expected1:difficulty2) -0.01
## cor(Intercept,phase1:difficulty1) -0.09
## cor(expected1,phase1:difficulty1) 0.06
## cor(phase1,phase1:difficulty1) 0.15
## cor(phase2,phase1:difficulty1) -0.11
## cor(difficulty1,phase1:difficulty1) 0.07
## cor(difficulty2,phase1:difficulty1) -0.06
## cor(expected1:phase1,phase1:difficulty1) 0.07
## cor(expected1:phase2,phase1:difficulty1) -0.04
## cor(expected1:difficulty1,phase1:difficulty1) 0.05
## cor(expected1:difficulty2,phase1:difficulty1) 0.01
## cor(Intercept,phase2:difficulty1) -0.05
## cor(expected1,phase2:difficulty1) 0.05
## cor(phase1,phase2:difficulty1) -0.11
## cor(phase2,phase2:difficulty1) 0.03
## cor(difficulty1,phase2:difficulty1) 0.09
## cor(difficulty2,phase2:difficulty1) 0.01
## cor(expected1:phase1,phase2:difficulty1) 0.01
## cor(expected1:phase2,phase2:difficulty1) 0.08
## cor(expected1:difficulty1,phase2:difficulty1) 0.03
## cor(expected1:difficulty2,phase2:difficulty1) -0.03
## cor(phase1:difficulty1,phase2:difficulty1) -0.02
## cor(Intercept,phase1:difficulty2) 0.06
## cor(expected1,phase1:difficulty2) -0.03
## cor(phase1,phase1:difficulty2) 0.04
## cor(phase2,phase1:difficulty2) -0.06
## cor(difficulty1,phase1:difficulty2) -0.05
## cor(difficulty2,phase1:difficulty2) 0.01
## cor(expected1:phase1,phase1:difficulty2) -0.03
## cor(expected1:phase2,phase1:difficulty2) -0.03
## cor(expected1:difficulty1,phase1:difficulty2) -0.04
## cor(expected1:difficulty2,phase1:difficulty2) -0.00
## cor(phase1:difficulty1,phase1:difficulty2) -0.05
## cor(phase2:difficulty1,phase1:difficulty2) -0.05
## cor(Intercept,phase2:difficulty2) 0.23
## cor(expected1,phase2:difficulty2) 0.13
## cor(phase1,phase2:difficulty2) -0.08
## cor(phase2,phase2:difficulty2) 0.04
## cor(difficulty1,phase2:difficulty2) 0.04
## cor(difficulty2,phase2:difficulty2) 0.02
## cor(expected1:phase1,phase2:difficulty2) 0.09
## cor(expected1:phase2,phase2:difficulty2) 0.19
## cor(expected1:difficulty1,phase2:difficulty2) -0.02
## cor(expected1:difficulty2,phase2:difficulty2) 0.01
## cor(phase1:difficulty1,phase2:difficulty2) -0.09
## cor(phase2:difficulty1,phase2:difficulty2) -0.04
## cor(phase1:difficulty2,phase2:difficulty2) -0.05
## cor(Intercept,expected1:phase1:difficulty1) 0.06
## cor(expected1,expected1:phase1:difficulty1) 0.07
## cor(phase1,expected1:phase1:difficulty1) 0.13
## cor(phase2,expected1:phase1:difficulty1) -0.01
## cor(difficulty1,expected1:phase1:difficulty1) -0.01
## cor(difficulty2,expected1:phase1:difficulty1) -0.01
## cor(expected1:phase1,expected1:phase1:difficulty1) 0.03
## cor(expected1:phase2,expected1:phase1:difficulty1) 0.02
## cor(expected1:difficulty1,expected1:phase1:difficulty1) 0.02
## cor(expected1:difficulty2,expected1:phase1:difficulty1) 0.04
## cor(phase1:difficulty1,expected1:phase1:difficulty1) -0.05
## cor(phase2:difficulty1,expected1:phase1:difficulty1) -0.04
## cor(phase1:difficulty2,expected1:phase1:difficulty1) -0.02
## cor(phase2:difficulty2,expected1:phase1:difficulty1) -0.00
## cor(Intercept,expected1:phase2:difficulty1) 0.04
## cor(expected1,expected1:phase2:difficulty1) 0.10
## cor(phase1,expected1:phase2:difficulty1) 0.06
## cor(phase2,expected1:phase2:difficulty1) 0.05
## cor(difficulty1,expected1:phase2:difficulty1) -0.04
## cor(difficulty2,expected1:phase2:difficulty1) -0.03
## cor(expected1:phase1,expected1:phase2:difficulty1) -0.01
## cor(expected1:phase2,expected1:phase2:difficulty1) -0.01
## cor(expected1:difficulty1,expected1:phase2:difficulty1) -0.02
## cor(expected1:difficulty2,expected1:phase2:difficulty1) 0.01
## cor(phase1:difficulty1,expected1:phase2:difficulty1) -0.02
## cor(phase2:difficulty1,expected1:phase2:difficulty1) -0.05
## cor(phase1:difficulty2,expected1:phase2:difficulty1) -0.01
## cor(phase2:difficulty2,expected1:phase2:difficulty1) -0.01
## cor(expected1:phase1:difficulty1,expected1:phase2:difficulty1) -0.00
## cor(Intercept,expected1:phase1:difficulty2) 0.02
## cor(expected1,expected1:phase1:difficulty2) -0.04
## cor(phase1,expected1:phase1:difficulty2) 0.02
## cor(phase2,expected1:phase1:difficulty2) -0.01
## cor(difficulty1,expected1:phase1:difficulty2) -0.02
## cor(difficulty2,expected1:phase1:difficulty2) 0.03
## cor(expected1:phase1,expected1:phase1:difficulty2) -0.03
## cor(expected1:phase2,expected1:phase1:difficulty2) -0.03
## cor(expected1:difficulty1,expected1:phase1:difficulty2) -0.02
## cor(expected1:difficulty2,expected1:phase1:difficulty2) 0.01
## cor(phase1:difficulty1,expected1:phase1:difficulty2) -0.04
## cor(phase2:difficulty1,expected1:phase1:difficulty2) -0.03
## cor(phase1:difficulty2,expected1:phase1:difficulty2) -0.01
## cor(phase2:difficulty2,expected1:phase1:difficulty2) -0.03
## cor(expected1:phase1:difficulty1,expected1:phase1:difficulty2) -0.03
## cor(expected1:phase2:difficulty1,expected1:phase1:difficulty2) -0.01
## cor(Intercept,expected1:phase2:difficulty2) 0.04
## cor(expected1,expected1:phase2:difficulty2) 0.07
## cor(phase1,expected1:phase2:difficulty2) 0.09
## cor(phase2,expected1:phase2:difficulty2) 0.06
## cor(difficulty1,expected1:phase2:difficulty2) -0.03
## cor(difficulty2,expected1:phase2:difficulty2) -0.04
## cor(expected1:phase1,expected1:phase2:difficulty2) 0.01
## cor(expected1:phase2,expected1:phase2:difficulty2) 0.01
## cor(expected1:difficulty1,expected1:phase2:difficulty2) -0.01
## cor(expected1:difficulty2,expected1:phase2:difficulty2) -0.00
## cor(phase1:difficulty1,expected1:phase2:difficulty2) -0.00
## cor(phase2:difficulty1,expected1:phase2:difficulty2) -0.04
## cor(phase1:difficulty2,expected1:phase2:difficulty2) -0.03
## cor(phase2:difficulty2,expected1:phase2:difficulty2) -0.07
## cor(expected1:phase1:difficulty1,expected1:phase2:difficulty2) 0.01
## cor(expected1:phase2:difficulty1,expected1:phase2:difficulty2) -0.01
## cor(expected1:phase1:difficulty2,expected1:phase2:difficulty2) -0.03
## Est.Error
## sd(Intercept) 0.02
## sd(expected1) 0.00
## sd(phase1) 0.01
## sd(phase2) 0.00
## sd(difficulty1) 0.00
## sd(difficulty2) 0.00
## sd(expected1:phase1) 0.01
## sd(expected1:phase2) 0.01
## sd(expected1:difficulty1) 0.00
## sd(expected1:difficulty2) 0.00
## sd(phase1:difficulty1) 0.01
## sd(phase2:difficulty1) 0.00
## sd(phase1:difficulty2) 0.00
## sd(phase2:difficulty2) 0.01
## sd(expected1:phase1:difficulty1) 0.00
## sd(expected1:phase2:difficulty1) 0.00
## sd(expected1:phase1:difficulty2) 0.00
## sd(expected1:phase2:difficulty2) 0.00
## cor(Intercept,expected1) 0.14
## cor(Intercept,phase1) 0.12
## cor(expected1,phase1) 0.15
## cor(Intercept,phase2) 0.14
## cor(expected1,phase2) 0.18
## cor(phase1,phase2) 0.15
## cor(Intercept,difficulty1) 0.17
## cor(expected1,difficulty1) 0.19
## cor(phase1,difficulty1) 0.17
## cor(phase2,difficulty1) 0.19
## cor(Intercept,difficulty2) 0.19
## cor(expected1,difficulty2) 0.20
## cor(phase1,difficulty2) 0.19
## cor(phase2,difficulty2) 0.20
## cor(difficulty1,difficulty2) 0.21
## cor(Intercept,expected1:phase1) 0.17
## cor(expected1,expected1:phase1) 0.19
## cor(phase1,expected1:phase1) 0.18
## cor(phase2,expected1:phase1) 0.19
## cor(difficulty1,expected1:phase1) 0.20
## cor(difficulty2,expected1:phase1) 0.21
## cor(Intercept,expected1:phase2) 0.17
## cor(expected1,expected1:phase2) 0.19
## cor(phase1,expected1:phase2) 0.17
## cor(phase2,expected1:phase2) 0.19
## cor(difficulty1,expected1:phase2) 0.20
## cor(difficulty2,expected1:phase2) 0.21
## cor(expected1:phase1,expected1:phase2) 0.20
## cor(Intercept,expected1:difficulty1) 0.21
## cor(expected1,expected1:difficulty1) 0.22
## cor(phase1,expected1:difficulty1) 0.21
## cor(phase2,expected1:difficulty1) 0.21
## cor(difficulty1,expected1:difficulty1) 0.22
## cor(difficulty2,expected1:difficulty1) 0.21
## cor(expected1:phase1,expected1:difficulty1) 0.22
## cor(expected1:phase2,expected1:difficulty1) 0.21
## cor(Intercept,expected1:difficulty2) 0.19
## cor(expected1,expected1:difficulty2) 0.21
## cor(phase1,expected1:difficulty2) 0.20
## cor(phase2,expected1:difficulty2) 0.21
## cor(difficulty1,expected1:difficulty2) 0.21
## cor(difficulty2,expected1:difficulty2) 0.22
## cor(expected1:phase1,expected1:difficulty2) 0.21
## cor(expected1:phase2,expected1:difficulty2) 0.21
## cor(expected1:difficulty1,expected1:difficulty2) 0.22
## cor(Intercept,phase1:difficulty1) 0.20
## cor(expected1,phase1:difficulty1) 0.21
## cor(phase1,phase1:difficulty1) 0.21
## cor(phase2,phase1:difficulty1) 0.21
## cor(difficulty1,phase1:difficulty1) 0.21
## cor(difficulty2,phase1:difficulty1) 0.21
## cor(expected1:phase1,phase1:difficulty1) 0.21
## cor(expected1:phase2,phase1:difficulty1) 0.21
## cor(expected1:difficulty1,phase1:difficulty1) 0.22
## cor(expected1:difficulty2,phase1:difficulty1) 0.22
## cor(Intercept,phase2:difficulty1) 0.20
## cor(expected1,phase2:difficulty1) 0.21
## cor(phase1,phase2:difficulty1) 0.20
## cor(phase2,phase2:difficulty1) 0.21
## cor(difficulty1,phase2:difficulty1) 0.21
## cor(difficulty2,phase2:difficulty1) 0.21
## cor(expected1:phase1,phase2:difficulty1) 0.21
## cor(expected1:phase2,phase2:difficulty1) 0.21
## cor(expected1:difficulty1,phase2:difficulty1) 0.22
## cor(expected1:difficulty2,phase2:difficulty1) 0.22
## cor(phase1:difficulty1,phase2:difficulty1) 0.22
## cor(Intercept,phase1:difficulty2) 0.21
## cor(expected1,phase1:difficulty2) 0.21
## cor(phase1,phase1:difficulty2) 0.21
## cor(phase2,phase1:difficulty2) 0.21
## cor(difficulty1,phase1:difficulty2) 0.22
## cor(difficulty2,phase1:difficulty2) 0.22
## cor(expected1:phase1,phase1:difficulty2) 0.22
## cor(expected1:phase2,phase1:difficulty2) 0.22
## cor(expected1:difficulty1,phase1:difficulty2) 0.22
## cor(expected1:difficulty2,phase1:difficulty2) 0.22
## cor(phase1:difficulty1,phase1:difficulty2) 0.22
## cor(phase2:difficulty1,phase1:difficulty2) 0.22
## cor(Intercept,phase2:difficulty2) 0.17
## cor(expected1,phase2:difficulty2) 0.19
## cor(phase1,phase2:difficulty2) 0.18
## cor(phase2,phase2:difficulty2) 0.19
## cor(difficulty1,phase2:difficulty2) 0.20
## cor(difficulty2,phase2:difficulty2) 0.21
## cor(expected1:phase1,phase2:difficulty2) 0.20
## cor(expected1:phase2,phase2:difficulty2) 0.20
## cor(expected1:difficulty1,phase2:difficulty2) 0.22
## cor(expected1:difficulty2,phase2:difficulty2) 0.21
## cor(phase1:difficulty1,phase2:difficulty2) 0.22
## cor(phase2:difficulty1,phase2:difficulty2) 0.21
## cor(phase1:difficulty2,phase2:difficulty2) 0.22
## cor(Intercept,expected1:phase1:difficulty1) 0.21
## cor(expected1,expected1:phase1:difficulty1) 0.22
## cor(phase1,expected1:phase1:difficulty1) 0.22
## cor(phase2,expected1:phase1:difficulty1) 0.22
## cor(difficulty1,expected1:phase1:difficulty1) 0.21
## cor(difficulty2,expected1:phase1:difficulty1) 0.22
## cor(expected1:phase1,expected1:phase1:difficulty1) 0.22
## cor(expected1:phase2,expected1:phase1:difficulty1) 0.21
## cor(expected1:difficulty1,expected1:phase1:difficulty1) 0.22
## cor(expected1:difficulty2,expected1:phase1:difficulty1) 0.22
## cor(phase1:difficulty1,expected1:phase1:difficulty1) 0.22
## cor(phase2:difficulty1,expected1:phase1:difficulty1) 0.22
## cor(phase1:difficulty2,expected1:phase1:difficulty1) 0.22
## cor(phase2:difficulty2,expected1:phase1:difficulty1) 0.22
## cor(Intercept,expected1:phase2:difficulty1) 0.21
## cor(expected1,expected1:phase2:difficulty1) 0.22
## cor(phase1,expected1:phase2:difficulty1) 0.21
## cor(phase2,expected1:phase2:difficulty1) 0.21
## cor(difficulty1,expected1:phase2:difficulty1) 0.22
## cor(difficulty2,expected1:phase2:difficulty1) 0.22
## cor(expected1:phase1,expected1:phase2:difficulty1) 0.22
## cor(expected1:phase2,expected1:phase2:difficulty1) 0.22
## cor(expected1:difficulty1,expected1:phase2:difficulty1) 0.22
## cor(expected1:difficulty2,expected1:phase2:difficulty1) 0.22
## cor(phase1:difficulty1,expected1:phase2:difficulty1) 0.22
## cor(phase2:difficulty1,expected1:phase2:difficulty1) 0.22
## cor(phase1:difficulty2,expected1:phase2:difficulty1) 0.22
## cor(phase2:difficulty2,expected1:phase2:difficulty1) 0.22
## cor(expected1:phase1:difficulty1,expected1:phase2:difficulty1) 0.22
## cor(Intercept,expected1:phase1:difficulty2) 0.21
## cor(expected1,expected1:phase1:difficulty2) 0.21
## cor(phase1,expected1:phase1:difficulty2) 0.21
## cor(phase2,expected1:phase1:difficulty2) 0.21
## cor(difficulty1,expected1:phase1:difficulty2) 0.21
## cor(difficulty2,expected1:phase1:difficulty2) 0.22
## cor(expected1:phase1,expected1:phase1:difficulty2) 0.21
## cor(expected1:phase2,expected1:phase1:difficulty2) 0.22
## cor(expected1:difficulty1,expected1:phase1:difficulty2) 0.22
## cor(expected1:difficulty2,expected1:phase1:difficulty2) 0.22
## cor(phase1:difficulty1,expected1:phase1:difficulty2) 0.22
## cor(phase2:difficulty1,expected1:phase1:difficulty2) 0.22
## cor(phase1:difficulty2,expected1:phase1:difficulty2) 0.22
## cor(phase2:difficulty2,expected1:phase1:difficulty2) 0.22
## cor(expected1:phase1:difficulty1,expected1:phase1:difficulty2) 0.22
## cor(expected1:phase2:difficulty1,expected1:phase1:difficulty2) 0.21
## cor(Intercept,expected1:phase2:difficulty2) 0.21
## cor(expected1,expected1:phase2:difficulty2) 0.21
## cor(phase1,expected1:phase2:difficulty2) 0.21
## cor(phase2,expected1:phase2:difficulty2) 0.21
## cor(difficulty1,expected1:phase2:difficulty2) 0.21
## cor(difficulty2,expected1:phase2:difficulty2) 0.22
## cor(expected1:phase1,expected1:phase2:difficulty2) 0.21
## cor(expected1:phase2,expected1:phase2:difficulty2) 0.21
## cor(expected1:difficulty1,expected1:phase2:difficulty2) 0.22
## cor(expected1:difficulty2,expected1:phase2:difficulty2) 0.22
## cor(phase1:difficulty1,expected1:phase2:difficulty2) 0.22
## cor(phase2:difficulty1,expected1:phase2:difficulty2) 0.22
## cor(phase1:difficulty2,expected1:phase2:difficulty2) 0.22
## cor(phase2:difficulty2,expected1:phase2:difficulty2) 0.22
## cor(expected1:phase1:difficulty1,expected1:phase2:difficulty2) 0.22
## cor(expected1:phase2:difficulty1,expected1:phase2:difficulty2) 0.22
## cor(expected1:phase1:difficulty2,expected1:phase2:difficulty2) 0.22
## l-95% CI
## sd(Intercept) 0.16
## sd(expected1) 0.01
## sd(phase1) 0.05
## sd(phase2) 0.01
## sd(difficulty1) 0.00
## sd(difficulty2) 0.00
## sd(expected1:phase1) 0.00
## sd(expected1:phase2) 0.00
## sd(expected1:difficulty1) 0.00
## sd(expected1:difficulty2) 0.00
## sd(phase1:difficulty1) 0.00
## sd(phase2:difficulty1) 0.00
## sd(phase1:difficulty2) 0.00
## sd(phase2:difficulty2) 0.00
## sd(expected1:phase1:difficulty1) 0.00
## sd(expected1:phase2:difficulty1) 0.00
## sd(expected1:phase1:difficulty2) 0.00
## sd(expected1:phase2:difficulty2) 0.00
## cor(Intercept,expected1) -0.06
## cor(Intercept,phase1) -0.15
## cor(expected1,phase1) -0.09
## cor(Intercept,phase2) -0.19
## cor(expected1,phase2) -0.31
## cor(phase1,phase2) -0.52
## cor(Intercept,difficulty1) -0.55
## cor(expected1,difficulty1) -0.44
## cor(phase1,difficulty1) -0.50
## cor(phase2,difficulty1) -0.44
## cor(Intercept,difficulty2) -0.40
## cor(expected1,difficulty2) -0.49
## cor(phase1,difficulty2) -0.52
## cor(phase2,difficulty2) -0.31
## cor(difficulty1,difficulty2) -0.38
## cor(Intercept,expected1:phase1) -0.30
## cor(expected1,expected1:phase1) -0.22
## cor(phase1,expected1:phase1) -0.25
## cor(phase2,expected1:phase1) -0.53
## cor(difficulty1,expected1:phase1) -0.33
## cor(difficulty2,expected1:phase1) -0.42
## cor(Intercept,expected1:phase2) -0.14
## cor(expected1,expected1:phase2) -0.30
## cor(phase1,expected1:phase2) -0.38
## cor(phase2,expected1:phase2) -0.33
## cor(difficulty1,expected1:phase2) -0.35
## cor(difficulty2,expected1:phase2) -0.29
## cor(expected1:phase1,expected1:phase2) -0.38
## cor(Intercept,expected1:difficulty1) -0.47
## cor(expected1,expected1:difficulty1) -0.39
## cor(phase1,expected1:difficulty1) -0.40
## cor(phase2,expected1:difficulty1) -0.42
## cor(difficulty1,expected1:difficulty1) -0.43
## cor(difficulty2,expected1:difficulty1) -0.42
## cor(expected1:phase1,expected1:difficulty1) -0.39
## cor(expected1:phase2,expected1:difficulty1) -0.40
## cor(Intercept,expected1:difficulty2) -0.37
## cor(expected1,expected1:difficulty2) -0.36
## cor(phase1,expected1:difficulty2) -0.31
## cor(phase2,expected1:difficulty2) -0.34
## cor(difficulty1,expected1:difficulty2) -0.39
## cor(difficulty2,expected1:difficulty2) -0.42
## cor(expected1:phase1,expected1:difficulty2) -0.36
## cor(expected1:phase2,expected1:difficulty2) -0.34
## cor(expected1:difficulty1,expected1:difficulty2) -0.42
## cor(Intercept,phase1:difficulty1) -0.46
## cor(expected1,phase1:difficulty1) -0.36
## cor(phase1,phase1:difficulty1) -0.28
## cor(phase2,phase1:difficulty1) -0.51
## cor(difficulty1,phase1:difficulty1) -0.35
## cor(difficulty2,phase1:difficulty1) -0.48
## cor(expected1:phase1,phase1:difficulty1) -0.35
## cor(expected1:phase2,phase1:difficulty1) -0.44
## cor(expected1:difficulty1,phase1:difficulty1) -0.38
## cor(expected1:difficulty2,phase1:difficulty1) -0.41
## cor(Intercept,phase2:difficulty1) -0.42
## cor(expected1,phase2:difficulty1) -0.37
## cor(phase1,phase2:difficulty1) -0.49
## cor(phase2,phase2:difficulty1) -0.37
## cor(difficulty1,phase2:difficulty1) -0.34
## cor(difficulty2,phase2:difficulty1) -0.40
## cor(expected1:phase1,phase2:difficulty1) -0.40
## cor(expected1:phase2,phase2:difficulty1) -0.35
## cor(expected1:difficulty1,phase2:difficulty1) -0.40
## cor(expected1:difficulty2,phase2:difficulty1) -0.45
## cor(phase1:difficulty1,phase2:difficulty1) -0.44
## cor(Intercept,phase1:difficulty2) -0.36
## cor(expected1,phase1:difficulty2) -0.43
## cor(phase1,phase1:difficulty2) -0.37
## cor(phase2,phase1:difficulty2) -0.46
## cor(difficulty1,phase1:difficulty2) -0.46
## cor(difficulty2,phase1:difficulty2) -0.41
## cor(expected1:phase1,phase1:difficulty2) -0.45
## cor(expected1:phase2,phase1:difficulty2) -0.44
## cor(expected1:difficulty1,phase1:difficulty2) -0.46
## cor(expected1:difficulty2,phase1:difficulty2) -0.42
## cor(phase1:difficulty1,phase1:difficulty2) -0.48
## cor(phase2:difficulty1,phase1:difficulty2) -0.47
## cor(Intercept,phase2:difficulty2) -0.14
## cor(expected1,phase2:difficulty2) -0.26
## cor(phase1,phase2:difficulty2) -0.42
## cor(phase2,phase2:difficulty2) -0.33
## cor(difficulty1,phase2:difficulty2) -0.36
## cor(difficulty2,phase2:difficulty2) -0.38
## cor(expected1:phase1,phase2:difficulty2) -0.32
## cor(expected1:phase2,phase2:difficulty2) -0.22
## cor(expected1:difficulty1,phase2:difficulty2) -0.43
## cor(expected1:difficulty2,phase2:difficulty2) -0.41
## cor(phase1:difficulty1,phase2:difficulty2) -0.49
## cor(phase2:difficulty1,phase2:difficulty2) -0.44
## cor(phase1:difficulty2,phase2:difficulty2) -0.46
## cor(Intercept,expected1:phase1:difficulty1) -0.36
## cor(expected1,expected1:phase1:difficulty1) -0.37
## cor(phase1,expected1:phase1:difficulty1) -0.33
## cor(phase2,expected1:phase1:difficulty1) -0.43
## cor(difficulty1,expected1:phase1:difficulty1) -0.42
## cor(difficulty2,expected1:phase1:difficulty1) -0.43
## cor(expected1:phase1,expected1:phase1:difficulty1) -0.40
## cor(expected1:phase2,expected1:phase1:difficulty1) -0.39
## cor(expected1:difficulty1,expected1:phase1:difficulty1) -0.41
## cor(expected1:difficulty2,expected1:phase1:difficulty1) -0.40
## cor(phase1:difficulty1,expected1:phase1:difficulty1) -0.47
## cor(phase2:difficulty1,expected1:phase1:difficulty1) -0.46
## cor(phase1:difficulty2,expected1:phase1:difficulty1) -0.44
## cor(phase2:difficulty2,expected1:phase1:difficulty1) -0.41
## cor(Intercept,expected1:phase2:difficulty1) -0.37
## cor(expected1,expected1:phase2:difficulty1) -0.34
## cor(phase1,expected1:phase2:difficulty1) -0.37
## cor(phase2,expected1:phase2:difficulty1) -0.37
## cor(difficulty1,expected1:phase2:difficulty1) -0.45
## cor(difficulty2,expected1:phase2:difficulty1) -0.44
## cor(expected1:phase1,expected1:phase2:difficulty1) -0.44
## cor(expected1:phase2,expected1:phase2:difficulty1) -0.42
## cor(expected1:difficulty1,expected1:phase2:difficulty1) -0.43
## cor(expected1:difficulty2,expected1:phase2:difficulty1) -0.41
## cor(phase1:difficulty1,expected1:phase2:difficulty1) -0.44
## cor(phase2:difficulty1,expected1:phase2:difficulty1) -0.47
## cor(phase1:difficulty2,expected1:phase2:difficulty1) -0.44
## cor(phase2:difficulty2,expected1:phase2:difficulty1) -0.42
## cor(expected1:phase1:difficulty1,expected1:phase2:difficulty1) -0.42
## cor(Intercept,expected1:phase1:difficulty2) -0.38
## cor(expected1,expected1:phase1:difficulty2) -0.45
## cor(phase1,expected1:phase1:difficulty2) -0.39
## cor(phase2,expected1:phase1:difficulty2) -0.43
## cor(difficulty1,expected1:phase1:difficulty2) -0.43
## cor(difficulty2,expected1:phase1:difficulty2) -0.40
## cor(expected1:phase1,expected1:phase1:difficulty2) -0.44
## cor(expected1:phase2,expected1:phase1:difficulty2) -0.45
## cor(expected1:difficulty1,expected1:phase1:difficulty2) -0.43
## cor(expected1:difficulty2,expected1:phase1:difficulty2) -0.41
## cor(phase1:difficulty1,expected1:phase1:difficulty2) -0.46
## cor(phase2:difficulty1,expected1:phase1:difficulty2) -0.45
## cor(phase1:difficulty2,expected1:phase1:difficulty2) -0.44
## cor(phase2:difficulty2,expected1:phase1:difficulty2) -0.45
## cor(expected1:phase1:difficulty1,expected1:phase1:difficulty2) -0.46
## cor(expected1:phase2:difficulty1,expected1:phase1:difficulty2) -0.43
## cor(Intercept,expected1:phase2:difficulty2) -0.37
## cor(expected1,expected1:phase2:difficulty2) -0.35
## cor(phase1,expected1:phase2:difficulty2) -0.33
## cor(phase2,expected1:phase2:difficulty2) -0.36
## cor(difficulty1,expected1:phase2:difficulty2) -0.44
## cor(difficulty2,expected1:phase2:difficulty2) -0.46
## cor(expected1:phase1,expected1:phase2:difficulty2) -0.40
## cor(expected1:phase2,expected1:phase2:difficulty2) -0.42
## cor(expected1:difficulty1,expected1:phase2:difficulty2) -0.42
## cor(expected1:difficulty2,expected1:phase2:difficulty2) -0.42
## cor(phase1:difficulty1,expected1:phase2:difficulty2) -0.42
## cor(phase2:difficulty1,expected1:phase2:difficulty2) -0.46
## cor(phase1:difficulty2,expected1:phase2:difficulty2) -0.45
## cor(phase2:difficulty2,expected1:phase2:difficulty2) -0.49
## cor(expected1:phase1:difficulty1,expected1:phase2:difficulty2) -0.41
## cor(expected1:phase2:difficulty1,expected1:phase2:difficulty2) -0.42
## cor(expected1:phase1:difficulty2,expected1:phase2:difficulty2) -0.46
## u-95% CI Rhat
## sd(Intercept) 0.23 1.00
## sd(expected1) 0.03 1.00
## sd(phase1) 0.07 1.00
## sd(phase2) 0.03 1.00
## sd(difficulty1) 0.02 1.00
## sd(difficulty2) 0.02 1.00
## sd(expected1:phase1) 0.03 1.00
## sd(expected1:phase2) 0.02 1.00
## sd(expected1:difficulty1) 0.01 1.00
## sd(expected1:difficulty2) 0.02 1.00
## sd(phase1:difficulty1) 0.02 1.00
## sd(phase2:difficulty1) 0.02 1.00
## sd(phase1:difficulty2) 0.02 1.00
## sd(phase2:difficulty2) 0.02 1.00
## sd(expected1:phase1:difficulty1) 0.02 1.00
## sd(expected1:phase2:difficulty1) 0.01 1.00
## sd(expected1:phase1:difficulty2) 0.02 1.00
## sd(expected1:phase2:difficulty2) 0.02 1.00
## cor(Intercept,expected1) 0.50 1.00
## cor(Intercept,phase1) 0.31 1.00
## cor(expected1,phase1) 0.51 1.00
## cor(Intercept,phase2) 0.36 1.00
## cor(expected1,phase2) 0.38 1.00
## cor(phase1,phase2) 0.05 1.00
## cor(Intercept,difficulty1) 0.10 1.00
## cor(expected1,difficulty1) 0.30 1.00
## cor(phase1,difficulty1) 0.18 1.00
## cor(phase2,difficulty1) 0.30 1.00
## cor(Intercept,difficulty2) 0.34 1.00
## cor(expected1,difficulty2) 0.30 1.00
## cor(phase1,difficulty2) 0.23 1.00
## cor(phase2,difficulty2) 0.47 1.00
## cor(difficulty1,difficulty2) 0.43 1.00
## cor(Intercept,expected1:phase1) 0.36 1.00
## cor(expected1,expected1:phase1) 0.53 1.00
## cor(phase1,expected1:phase1) 0.44 1.00
## cor(phase2,expected1:phase1) 0.21 1.00
## cor(difficulty1,expected1:phase1) 0.44 1.00
## cor(difficulty2,expected1:phase1) 0.38 1.00
## cor(Intercept,expected1:phase2) 0.51 1.00
## cor(expected1,expected1:phase2) 0.45 1.00
## cor(phase1,expected1:phase2) 0.30 1.00
## cor(phase2,expected1:phase2) 0.43 1.00
## cor(difficulty1,expected1:phase2) 0.43 1.00
## cor(difficulty2,expected1:phase2) 0.52 1.00
## cor(expected1:phase1,expected1:phase2) 0.41 1.00
## cor(Intercept,expected1:difficulty1) 0.35 1.00
## cor(expected1,expected1:difficulty1) 0.45 1.00
## cor(phase1,expected1:difficulty1) 0.43 1.00
## cor(phase2,expected1:difficulty1) 0.42 1.00
## cor(difficulty1,expected1:difficulty1) 0.41 1.00
## cor(difficulty2,expected1:difficulty1) 0.41 1.00
## cor(expected1:phase1,expected1:difficulty1) 0.45 1.00
## cor(expected1:phase2,expected1:difficulty1) 0.43 1.00
## cor(Intercept,expected1:difficulty2) 0.38 1.00
## cor(expected1,expected1:difficulty2) 0.45 1.00
## cor(phase1,expected1:difficulty2) 0.47 1.00
## cor(phase2,expected1:difficulty2) 0.46 1.00
## cor(difficulty1,expected1:difficulty2) 0.42 1.00
## cor(difficulty2,expected1:difficulty2) 0.41 1.00
## cor(expected1:phase1,expected1:difficulty2) 0.46 1.00
## cor(expected1:phase2,expected1:difficulty2) 0.48 1.00
## cor(expected1:difficulty1,expected1:difficulty2) 0.42 1.00
## cor(Intercept,phase1:difficulty1) 0.31 1.00
## cor(expected1,phase1:difficulty1) 0.45 1.00
## cor(phase1,phase1:difficulty1) 0.53 1.00
## cor(phase2,phase1:difficulty1) 0.31 1.00
## cor(difficulty1,phase1:difficulty1) 0.47 1.00
## cor(difficulty2,phase1:difficulty1) 0.36 1.00
## cor(expected1:phase1,phase1:difficulty1) 0.47 1.00
## cor(expected1:phase2,phase1:difficulty1) 0.38 1.00
## cor(expected1:difficulty1,phase1:difficulty1) 0.46 1.00
## cor(expected1:difficulty2,phase1:difficulty1) 0.42 1.00
## cor(Intercept,phase2:difficulty1) 0.35 1.00
## cor(expected1,phase2:difficulty1) 0.45 1.00
## cor(phase1,phase2:difficulty1) 0.31 1.00
## cor(phase2,phase2:difficulty1) 0.42 1.00
## cor(difficulty1,phase2:difficulty1) 0.49 1.00
## cor(difficulty2,phase2:difficulty1) 0.42 1.00
## cor(expected1:phase1,phase2:difficulty1) 0.43 1.00
## cor(expected1:phase2,phase2:difficulty1) 0.48 1.00
## cor(expected1:difficulty1,phase2:difficulty1) 0.44 1.00
## cor(expected1:difficulty2,phase2:difficulty1) 0.40 1.00
## cor(phase1:difficulty1,phase2:difficulty1) 0.40 1.00
## cor(Intercept,phase1:difficulty2) 0.45 1.00
## cor(expected1,phase1:difficulty2) 0.39 1.00
## cor(phase1,phase1:difficulty2) 0.44 1.00
## cor(phase2,phase1:difficulty2) 0.36 1.00
## cor(difficulty1,phase1:difficulty2) 0.39 1.00
## cor(difficulty2,phase1:difficulty2) 0.43 1.00
## cor(expected1:phase1,phase1:difficulty2) 0.39 1.00
## cor(expected1:phase2,phase1:difficulty2) 0.38 1.00
## cor(expected1:difficulty1,phase1:difficulty2) 0.40 1.00
## cor(expected1:difficulty2,phase1:difficulty2) 0.42 1.00
## cor(phase1:difficulty1,phase1:difficulty2) 0.38 1.00
## cor(phase2:difficulty1,phase1:difficulty2) 0.38 1.00
## cor(Intercept,phase2:difficulty2) 0.54 1.00
## cor(expected1,phase2:difficulty2) 0.50 1.00
## cor(phase1,phase2:difficulty2) 0.28 1.00
## cor(phase2,phase2:difficulty2) 0.42 1.00
## cor(difficulty1,phase2:difficulty2) 0.42 1.00
## cor(difficulty2,phase2:difficulty2) 0.42 1.00
## cor(expected1:phase1,phase2:difficulty2) 0.47 1.00
## cor(expected1:phase2,phase2:difficulty2) 0.56 1.00
## cor(expected1:difficulty1,phase2:difficulty2) 0.41 1.00
## cor(expected1:difficulty2,phase2:difficulty2) 0.41 1.00
## cor(phase1:difficulty1,phase2:difficulty2) 0.35 1.00
## cor(phase2:difficulty1,phase2:difficulty2) 0.39 1.00
## cor(phase1:difficulty2,phase2:difficulty2) 0.38 1.00
## cor(Intercept,expected1:phase1:difficulty1) 0.46 1.00
## cor(expected1,expected1:phase1:difficulty1) 0.47 1.00
## cor(phase1,expected1:phase1:difficulty1) 0.54 1.00
## cor(phase2,expected1:phase1:difficulty1) 0.41 1.00
## cor(difficulty1,expected1:phase1:difficulty1) 0.41 1.00
## cor(difficulty2,expected1:phase1:difficulty1) 0.40 1.00
## cor(expected1:phase1,expected1:phase1:difficulty1) 0.45 1.00
## cor(expected1:phase2,expected1:phase1:difficulty1) 0.43 1.00
## cor(expected1:difficulty1,expected1:phase1:difficulty1) 0.44 1.00
## cor(expected1:difficulty2,expected1:phase1:difficulty1) 0.46 1.00
## cor(phase1:difficulty1,expected1:phase1:difficulty1) 0.39 1.00
## cor(phase2:difficulty1,expected1:phase1:difficulty1) 0.39 1.00
## cor(phase1:difficulty2,expected1:phase1:difficulty1) 0.41 1.00
## cor(phase2:difficulty2,expected1:phase1:difficulty1) 0.42 1.00
## cor(Intercept,expected1:phase2:difficulty1) 0.44 1.00
## cor(expected1,expected1:phase2:difficulty1) 0.50 1.00
## cor(phase1,expected1:phase2:difficulty1) 0.46 1.00
## cor(phase2,expected1:phase2:difficulty1) 0.46 1.00
## cor(difficulty1,expected1:phase2:difficulty1) 0.38 1.00
## cor(difficulty2,expected1:phase2:difficulty1) 0.39 1.00
## cor(expected1:phase1,expected1:phase2:difficulty1) 0.41 1.00
## cor(expected1:phase2,expected1:phase2:difficulty1) 0.41 1.00
## cor(expected1:difficulty1,expected1:phase2:difficulty1) 0.41 1.00
## cor(expected1:difficulty2,expected1:phase2:difficulty1) 0.42 1.00
## cor(phase1:difficulty1,expected1:phase2:difficulty1) 0.40 1.00
## cor(phase2:difficulty1,expected1:phase2:difficulty1) 0.38 1.00
## cor(phase1:difficulty2,expected1:phase2:difficulty1) 0.42 1.00
## cor(phase2:difficulty2,expected1:phase2:difficulty1) 0.41 1.00
## cor(expected1:phase1:difficulty1,expected1:phase2:difficulty1) 0.42 1.00
## cor(Intercept,expected1:phase1:difficulty2) 0.42 1.00
## cor(expected1,expected1:phase1:difficulty2) 0.39 1.00
## cor(phase1,expected1:phase1:difficulty2) 0.42 1.00
## cor(phase2,expected1:phase1:difficulty2) 0.40 1.00
## cor(difficulty1,expected1:phase1:difficulty2) 0.39 1.00
## cor(difficulty2,expected1:phase1:difficulty2) 0.44 1.00
## cor(expected1:phase1,expected1:phase1:difficulty2) 0.39 1.00
## cor(expected1:phase2,expected1:phase1:difficulty2) 0.40 1.00
## cor(expected1:difficulty1,expected1:phase1:difficulty2) 0.41 1.00
## cor(expected1:difficulty2,expected1:phase1:difficulty2) 0.43 1.00
## cor(phase1:difficulty1,expected1:phase1:difficulty2) 0.39 1.00
## cor(phase2:difficulty1,expected1:phase1:difficulty2) 0.40 1.00
## cor(phase1:difficulty2,expected1:phase1:difficulty2) 0.41 1.00
## cor(phase2:difficulty2,expected1:phase1:difficulty2) 0.39 1.00
## cor(expected1:phase1:difficulty1,expected1:phase1:difficulty2) 0.41 1.00
## cor(expected1:phase2:difficulty1,expected1:phase1:difficulty2) 0.40 1.00
## cor(Intercept,expected1:phase2:difficulty2) 0.43 1.00
## cor(expected1,expected1:phase2:difficulty2) 0.48 1.00
## cor(phase1,expected1:phase2:difficulty2) 0.48 1.00
## cor(phase2,expected1:phase2:difficulty2) 0.46 1.00
## cor(difficulty1,expected1:phase2:difficulty2) 0.39 1.00
## cor(difficulty2,expected1:phase2:difficulty2) 0.38 1.00
## cor(expected1:phase1,expected1:phase2:difficulty2) 0.43 1.00
## cor(expected1:phase2,expected1:phase2:difficulty2) 0.42 1.00
## cor(expected1:difficulty1,expected1:phase2:difficulty2) 0.42 1.00
## cor(expected1:difficulty2,expected1:phase2:difficulty2) 0.42 1.00
## cor(phase1:difficulty1,expected1:phase2:difficulty2) 0.42 1.00
## cor(phase2:difficulty1,expected1:phase2:difficulty2) 0.39 1.00
## cor(phase1:difficulty2,expected1:phase2:difficulty2) 0.40 1.00
## cor(phase2:difficulty2,expected1:phase2:difficulty2) 0.37 1.00
## cor(expected1:phase1:difficulty1,expected1:phase2:difficulty2) 0.44 1.00
## cor(expected1:phase2:difficulty1,expected1:phase2:difficulty2) 0.42 1.00
## cor(expected1:phase1:difficulty2,expected1:phase2:difficulty2) 0.39 1.00
## Bulk_ESS
## sd(Intercept) 2366
## sd(expected1) 7066
## sd(phase1) 7239
## sd(phase2) 5858
## sd(difficulty1) 4134
## sd(difficulty2) 3095
## sd(expected1:phase1) 3029
## sd(expected1:phase2) 3068
## sd(expected1:difficulty1) 6251
## sd(expected1:difficulty2) 3722
## sd(phase1:difficulty1) 5189
## sd(phase2:difficulty1) 4898
## sd(phase1:difficulty2) 6330
## sd(phase2:difficulty2) 5160
## sd(expected1:phase1:difficulty1) 7427
## sd(expected1:phase2:difficulty1) 6694
## sd(expected1:phase1:difficulty2) 6824
## sd(expected1:phase2:difficulty2) 5243
## cor(Intercept,expected1) 14133
## cor(Intercept,phase1) 6683
## cor(expected1,phase1) 2663
## cor(Intercept,phase2) 12625
## cor(expected1,phase2) 6010
## cor(phase1,phase2) 10967
## cor(Intercept,difficulty1) 16432
## cor(expected1,difficulty1) 12974
## cor(phase1,difficulty1) 15222
## cor(phase2,difficulty1) 13511
## cor(Intercept,difficulty2) 18051
## cor(expected1,difficulty2) 13174
## cor(phase1,difficulty2) 14317
## cor(phase2,difficulty2) 15036
## cor(difficulty1,difficulty2) 15041
## cor(Intercept,expected1:phase1) 18603
## cor(expected1,expected1:phase1) 9527
## cor(phase1,expected1:phase1) 13295
## cor(phase2,expected1:phase1) 9736
## cor(difficulty1,expected1:phase1) 10090
## cor(difficulty2,expected1:phase1) 8434
## cor(Intercept,expected1:phase2) 15949
## cor(expected1,expected1:phase2) 11900
## cor(phase1,expected1:phase2) 14657
## cor(phase2,expected1:phase2) 11006
## cor(difficulty1,expected1:phase2) 9619
## cor(difficulty2,expected1:phase2) 6869
## cor(expected1:phase1,expected1:phase2) 9188
## cor(Intercept,expected1:difficulty1) 22552
## cor(expected1,expected1:difficulty1) 23500
## cor(phase1,expected1:difficulty1) 23198
## cor(phase2,expected1:difficulty1) 23837
## cor(difficulty1,expected1:difficulty1) 19936
## cor(difficulty2,expected1:difficulty1) 17646
## cor(expected1:phase1,expected1:difficulty1) 16596
## cor(expected1:phase2,expected1:difficulty1) 18736
## cor(Intercept,expected1:difficulty2) 20349
## cor(expected1,expected1:difficulty2) 15884
## cor(phase1,expected1:difficulty2) 18050
## cor(phase2,expected1:difficulty2) 16314
## cor(difficulty1,expected1:difficulty2) 15940
## cor(difficulty2,expected1:difficulty2) 15629
## cor(expected1:phase1,expected1:difficulty2) 13665
## cor(expected1:phase2,expected1:difficulty2) 13394
## cor(expected1:difficulty1,expected1:difficulty2) 12206
## cor(Intercept,phase1:difficulty1) 20678
## cor(expected1,phase1:difficulty1) 16601
## cor(phase1,phase1:difficulty1) 15476
## cor(phase2,phase1:difficulty1) 15607
## cor(difficulty1,phase1:difficulty1) 16586
## cor(difficulty2,phase1:difficulty1) 15840
## cor(expected1:phase1,phase1:difficulty1) 17038
## cor(expected1:phase2,phase1:difficulty1) 17604
## cor(expected1:difficulty1,phase1:difficulty1) 13039
## cor(expected1:difficulty2,phase1:difficulty1) 14447
## cor(Intercept,phase2:difficulty1) 21087
## cor(expected1,phase2:difficulty1) 18351
## cor(phase1,phase2:difficulty1) 19609
## cor(phase2,phase2:difficulty1) 19276
## cor(difficulty1,phase2:difficulty1) 14912
## cor(difficulty2,phase2:difficulty1) 16387
## cor(expected1:phase1,phase2:difficulty1) 16007
## cor(expected1:phase2,phase2:difficulty1) 14459
## cor(expected1:difficulty1,phase2:difficulty1) 13914
## cor(expected1:difficulty2,phase2:difficulty1) 14616
## cor(phase1:difficulty1,phase2:difficulty1) 14149
## cor(Intercept,phase1:difficulty2) 23002
## cor(expected1,phase1:difficulty2) 21002
## cor(phase1,phase1:difficulty2) 21507
## cor(phase2,phase1:difficulty2) 19294
## cor(difficulty1,phase1:difficulty2) 18449
## cor(difficulty2,phase1:difficulty2) 16976
## cor(expected1:phase1,phase1:difficulty2) 18590
## cor(expected1:phase2,phase1:difficulty2) 18237
## cor(expected1:difficulty1,phase1:difficulty2) 12947
## cor(expected1:difficulty2,phase1:difficulty2) 15489
## cor(phase1:difficulty1,phase1:difficulty2) 12117
## cor(phase2:difficulty1,phase1:difficulty2) 12498
## cor(Intercept,phase2:difficulty2) 18737
## cor(expected1,phase2:difficulty2) 11983
## cor(phase1,phase2:difficulty2) 19341
## cor(phase2,phase2:difficulty2) 14594
## cor(difficulty1,phase2:difficulty2) 13320
## cor(difficulty2,phase2:difficulty2) 11807
## cor(expected1:phase1,phase2:difficulty2) 12466
## cor(expected1:phase2,phase2:difficulty2) 10284
## cor(expected1:difficulty1,phase2:difficulty2) 11650
## cor(expected1:difficulty2,phase2:difficulty2) 12079
## cor(phase1:difficulty1,phase2:difficulty2) 11518
## cor(phase2:difficulty1,phase2:difficulty2) 11069
## cor(phase1:difficulty2,phase2:difficulty2) 11275
## cor(Intercept,expected1:phase1:difficulty1) 23070
## cor(expected1,expected1:phase1:difficulty1) 19927
## cor(phase1,expected1:phase1:difficulty1) 19826
## cor(phase2,expected1:phase1:difficulty1) 23301
## cor(difficulty1,expected1:phase1:difficulty1) 19693
## cor(difficulty2,expected1:phase1:difficulty1) 18241
## cor(expected1:phase1,expected1:phase1:difficulty1) 20540
## cor(expected1:phase2,expected1:phase1:difficulty1) 18918
## cor(expected1:difficulty1,expected1:phase1:difficulty1) 14901
## cor(expected1:difficulty2,expected1:phase1:difficulty1) 15363
## cor(phase1:difficulty1,expected1:phase1:difficulty1) 13534
## cor(phase2:difficulty1,expected1:phase1:difficulty1) 12602
## cor(phase1:difficulty2,expected1:phase1:difficulty1) 12516
## cor(phase2:difficulty2,expected1:phase1:difficulty1) 15330
## cor(Intercept,expected1:phase2:difficulty1) 23415
## cor(expected1,expected1:phase2:difficulty1) 19618
## cor(phase1,expected1:phase2:difficulty1) 22784
## cor(phase2,expected1:phase2:difficulty1) 20412
## cor(difficulty1,expected1:phase2:difficulty1) 19912
## cor(difficulty2,expected1:phase2:difficulty1) 17690
## cor(expected1:phase1,expected1:phase2:difficulty1) 17898
## cor(expected1:phase2,expected1:phase2:difficulty1) 19036
## cor(expected1:difficulty1,expected1:phase2:difficulty1) 14669
## cor(expected1:difficulty2,expected1:phase2:difficulty1) 16225
## cor(phase1:difficulty1,expected1:phase2:difficulty1) 14340
## cor(phase2:difficulty1,expected1:phase2:difficulty1) 12139
## cor(phase1:difficulty2,expected1:phase2:difficulty1) 12162
## cor(phase2:difficulty2,expected1:phase2:difficulty1) 13742
## cor(expected1:phase1:difficulty1,expected1:phase2:difficulty1) 11303
## cor(Intercept,expected1:phase1:difficulty2) 24977
## cor(expected1,expected1:phase1:difficulty2) 21984
## cor(phase1,expected1:phase1:difficulty2) 25112
## cor(phase2,expected1:phase1:difficulty2) 23846
## cor(difficulty1,expected1:phase1:difficulty2) 18818
## cor(difficulty2,expected1:phase1:difficulty2) 16288
## cor(expected1:phase1,expected1:phase1:difficulty2) 18312
## cor(expected1:phase2,expected1:phase1:difficulty2) 18233
## cor(expected1:difficulty1,expected1:phase1:difficulty2) 13677
## cor(expected1:difficulty2,expected1:phase1:difficulty2) 14578
## cor(phase1:difficulty1,expected1:phase1:difficulty2) 13673
## cor(phase2:difficulty1,expected1:phase1:difficulty2) 13249
## cor(phase1:difficulty2,expected1:phase1:difficulty2) 12764
## cor(phase2:difficulty2,expected1:phase1:difficulty2) 15601
## cor(expected1:phase1:difficulty1,expected1:phase1:difficulty2) 11587
## cor(expected1:phase2:difficulty1,expected1:phase1:difficulty2) 10509
## cor(Intercept,expected1:phase2:difficulty2) 22721
## cor(expected1,expected1:phase2:difficulty2) 18641
## cor(phase1,expected1:phase2:difficulty2) 21970
## cor(phase2,expected1:phase2:difficulty2) 17872
## cor(difficulty1,expected1:phase2:difficulty2) 19510
## cor(difficulty2,expected1:phase2:difficulty2) 16470
## cor(expected1:phase1,expected1:phase2:difficulty2) 19280
## cor(expected1:phase2,expected1:phase2:difficulty2) 17766
## cor(expected1:difficulty1,expected1:phase2:difficulty2) 14471
## cor(expected1:difficulty2,expected1:phase2:difficulty2) 14774
## cor(phase1:difficulty1,expected1:phase2:difficulty2) 15075
## cor(phase2:difficulty1,expected1:phase2:difficulty2) 13598
## cor(phase1:difficulty2,expected1:phase2:difficulty2) 12355
## cor(phase2:difficulty2,expected1:phase2:difficulty2) 13156
## cor(expected1:phase1:difficulty1,expected1:phase2:difficulty2) 11294
## cor(expected1:phase2:difficulty1,expected1:phase2:difficulty2) 11446
## cor(expected1:phase1:difficulty2,expected1:phase2:difficulty2) 10626
## Tail_ESS
## sd(Intercept) 4846
## sd(expected1) 7666
## sd(phase1) 11628
## sd(phase2) 7273
## sd(difficulty1) 2776
## sd(difficulty2) 3998
## sd(expected1:phase1) 2182
## sd(expected1:phase2) 2836
## sd(expected1:difficulty1) 7980
## sd(expected1:difficulty2) 5876
## sd(phase1:difficulty1) 6095
## sd(phase2:difficulty1) 6781
## sd(phase1:difficulty2) 6790
## sd(phase2:difficulty2) 5173
## sd(expected1:phase1:difficulty1) 6748
## sd(expected1:phase2:difficulty1) 7628
## sd(expected1:phase1:difficulty2) 7751
## sd(expected1:phase2:difficulty2) 7346
## cor(Intercept,expected1) 13901
## cor(Intercept,phase1) 10052
## cor(expected1,phase1) 5453
## cor(Intercept,phase2) 13120
## cor(expected1,phase2) 8180
## cor(phase1,phase2) 13010
## cor(Intercept,difficulty1) 10045
## cor(expected1,difficulty1) 13322
## cor(phase1,difficulty1) 11806
## cor(phase2,difficulty1) 12642
## cor(Intercept,difficulty2) 12744
## cor(expected1,difficulty2) 13916
## cor(phase1,difficulty2) 11596
## cor(phase2,difficulty2) 13506
## cor(difficulty1,difficulty2) 13240
## cor(Intercept,expected1:phase1) 13281
## cor(expected1,expected1:phase1) 11700
## cor(phase1,expected1:phase1) 12199
## cor(phase2,expected1:phase1) 12181
## cor(difficulty1,expected1:phase1) 11920
## cor(difficulty2,expected1:phase1) 13198
## cor(Intercept,expected1:phase2) 12123
## cor(expected1,expected1:phase2) 13146
## cor(phase1,expected1:phase2) 13417
## cor(phase2,expected1:phase2) 12579
## cor(difficulty1,expected1:phase2) 11882
## cor(difficulty2,expected1:phase2) 11126
## cor(expected1:phase1,expected1:phase2) 13013
## cor(Intercept,expected1:difficulty1) 13046
## cor(expected1,expected1:difficulty1) 13311
## cor(phase1,expected1:difficulty1) 14136
## cor(phase2,expected1:difficulty1) 13834
## cor(difficulty1,expected1:difficulty1) 13714
## cor(difficulty2,expected1:difficulty1) 14953
## cor(expected1:phase1,expected1:difficulty1) 14001
## cor(expected1:phase2,expected1:difficulty1) 14017
## cor(Intercept,expected1:difficulty2) 13541
## cor(expected1,expected1:difficulty2) 13431
## cor(phase1,expected1:difficulty2) 13038
## cor(phase2,expected1:difficulty2) 13622
## cor(difficulty1,expected1:difficulty2) 12615
## cor(difficulty2,expected1:difficulty2) 13899
## cor(expected1:phase1,expected1:difficulty2) 14302
## cor(expected1:phase2,expected1:difficulty2) 14382
## cor(expected1:difficulty1,expected1:difficulty2) 14845
## cor(Intercept,phase1:difficulty1) 13282
## cor(expected1,phase1:difficulty1) 13648
## cor(phase1,phase1:difficulty1) 12386
## cor(phase2,phase1:difficulty1) 13542
## cor(difficulty1,phase1:difficulty1) 13964
## cor(difficulty2,phase1:difficulty1) 14759
## cor(expected1:phase1,phase1:difficulty1) 14109
## cor(expected1:phase2,phase1:difficulty1) 13794
## cor(expected1:difficulty1,phase1:difficulty1) 14056
## cor(expected1:difficulty2,phase1:difficulty1) 13721
## cor(Intercept,phase2:difficulty1) 13027
## cor(expected1,phase2:difficulty1) 13619
## cor(phase1,phase2:difficulty1) 13105
## cor(phase2,phase2:difficulty1) 13993
## cor(difficulty1,phase2:difficulty1) 13596
## cor(difficulty2,phase2:difficulty1) 14732
## cor(expected1:phase1,phase2:difficulty1) 14957
## cor(expected1:phase2,phase2:difficulty1) 14145
## cor(expected1:difficulty1,phase2:difficulty1) 13741
## cor(expected1:difficulty2,phase2:difficulty1) 15255
## cor(phase1:difficulty1,phase2:difficulty1) 15542
## cor(Intercept,phase1:difficulty2) 13489
## cor(expected1,phase1:difficulty2) 13678
## cor(phase1,phase1:difficulty2) 12851
## cor(phase2,phase1:difficulty2) 13674
## cor(difficulty1,phase1:difficulty2) 13341
## cor(difficulty2,phase1:difficulty2) 13816
## cor(expected1:phase1,phase1:difficulty2) 13784
## cor(expected1:phase2,phase1:difficulty2) 14706
## cor(expected1:difficulty1,phase1:difficulty2) 13261
## cor(expected1:difficulty2,phase1:difficulty2) 14463
## cor(phase1:difficulty1,phase1:difficulty2) 14085
## cor(phase2:difficulty1,phase1:difficulty2) 14235
## cor(Intercept,phase2:difficulty2) 12832
## cor(expected1,phase2:difficulty2) 13495
## cor(phase1,phase2:difficulty2) 13723
## cor(phase2,phase2:difficulty2) 13388
## cor(difficulty1,phase2:difficulty2) 13522
## cor(difficulty2,phase2:difficulty2) 13222
## cor(expected1:phase1,phase2:difficulty2) 13607
## cor(expected1:phase2,phase2:difficulty2) 12021
## cor(expected1:difficulty1,phase2:difficulty2) 14430
## cor(expected1:difficulty2,phase2:difficulty2) 13639
## cor(phase1:difficulty1,phase2:difficulty2) 13764
## cor(phase2:difficulty1,phase2:difficulty2) 13464
## cor(phase1:difficulty2,phase2:difficulty2) 15136
## cor(Intercept,expected1:phase1:difficulty1) 13700
## cor(expected1,expected1:phase1:difficulty1) 13061
## cor(phase1,expected1:phase1:difficulty1) 13947
## cor(phase2,expected1:phase1:difficulty1) 13971
## cor(difficulty1,expected1:phase1:difficulty1) 14156
## cor(difficulty2,expected1:phase1:difficulty1) 13666
## cor(expected1:phase1,expected1:phase1:difficulty1) 13708
## cor(expected1:phase2,expected1:phase1:difficulty1) 14229
## cor(expected1:difficulty1,expected1:phase1:difficulty1) 14680
## cor(expected1:difficulty2,expected1:phase1:difficulty1) 14001
## cor(phase1:difficulty1,expected1:phase1:difficulty1) 13316
## cor(phase2:difficulty1,expected1:phase1:difficulty1) 14134
## cor(phase1:difficulty2,expected1:phase1:difficulty1) 14264
## cor(phase2:difficulty2,expected1:phase1:difficulty1) 15749
## cor(Intercept,expected1:phase2:difficulty1) 13120
## cor(expected1,expected1:phase2:difficulty1) 13801
## cor(phase1,expected1:phase2:difficulty1) 13974
## cor(phase2,expected1:phase2:difficulty1) 14136
## cor(difficulty1,expected1:phase2:difficulty1) 14253
## cor(difficulty2,expected1:phase2:difficulty1) 14454
## cor(expected1:phase1,expected1:phase2:difficulty1) 13891
## cor(expected1:phase2,expected1:phase2:difficulty1) 14646
## cor(expected1:difficulty1,expected1:phase2:difficulty1) 14325
## cor(expected1:difficulty2,expected1:phase2:difficulty1) 15700
## cor(phase1:difficulty1,expected1:phase2:difficulty1) 14899
## cor(phase2:difficulty1,expected1:phase2:difficulty1) 14768
## cor(phase1:difficulty2,expected1:phase2:difficulty1) 13772
## cor(phase2:difficulty2,expected1:phase2:difficulty1) 14138
## cor(expected1:phase1:difficulty1,expected1:phase2:difficulty1) 14170
## cor(Intercept,expected1:phase1:difficulty2) 13094
## cor(expected1,expected1:phase1:difficulty2) 13509
## cor(phase1,expected1:phase1:difficulty2) 12429
## cor(phase2,expected1:phase1:difficulty2) 14026
## cor(difficulty1,expected1:phase1:difficulty2) 14670
## cor(difficulty2,expected1:phase1:difficulty2) 14558
## cor(expected1:phase1,expected1:phase1:difficulty2) 14823
## cor(expected1:phase2,expected1:phase1:difficulty2) 14127
## cor(expected1:difficulty1,expected1:phase1:difficulty2) 14486
## cor(expected1:difficulty2,expected1:phase1:difficulty2) 14711
## cor(phase1:difficulty1,expected1:phase1:difficulty2) 14330
## cor(phase2:difficulty1,expected1:phase1:difficulty2) 14702
## cor(phase1:difficulty2,expected1:phase1:difficulty2) 15253
## cor(phase2:difficulty2,expected1:phase1:difficulty2) 15659
## cor(expected1:phase1:difficulty1,expected1:phase1:difficulty2) 14662
## cor(expected1:phase2:difficulty1,expected1:phase1:difficulty2) 13875
## cor(Intercept,expected1:phase2:difficulty2) 13616
## cor(expected1,expected1:phase2:difficulty2) 14173
## cor(phase1,expected1:phase2:difficulty2) 13712
## cor(phase2,expected1:phase2:difficulty2) 13197
## cor(difficulty1,expected1:phase2:difficulty2) 14104
## cor(difficulty2,expected1:phase2:difficulty2) 14255
## cor(expected1:phase1,expected1:phase2:difficulty2) 14359
## cor(expected1:phase2,expected1:phase2:difficulty2) 14830
## cor(expected1:difficulty1,expected1:phase2:difficulty2) 15253
## cor(expected1:difficulty2,expected1:phase2:difficulty2) 14510
## cor(phase1:difficulty1,expected1:phase2:difficulty2) 14737
## cor(phase2:difficulty1,expected1:phase2:difficulty2) 14498
## cor(phase1:difficulty2,expected1:phase2:difficulty2) 13653
## cor(phase2:difficulty2,expected1:phase2:difficulty2) 14074
## cor(expected1:phase1:difficulty1,expected1:phase2:difficulty2) 14532
## cor(expected1:phase2:difficulty1,expected1:phase2:difficulty2) 14123
## cor(expected1:phase1:difficulty2,expected1:phase2:difficulty2) 14521
##
## ~trl (Number of levels: 288)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## sd(Intercept) 0.07 0.00 0.06 0.08 1.00 6048
## sd(diagnosis1) 0.01 0.00 0.00 0.02 1.00 3418
## sd(diagnosis2) 0.00 0.00 0.00 0.01 1.00 4559
## cor(Intercept,diagnosis1) 0.49 0.23 -0.03 0.87 1.00 10682
## cor(Intercept,diagnosis2) -0.03 0.37 -0.72 0.69 1.00 18822
## cor(diagnosis1,diagnosis2) -0.23 0.42 -0.88 0.65 1.00 8090
## Tail_ESS
## sd(Intercept) 10131
## sd(diagnosis1) 3325
## sd(diagnosis2) 7185
## cor(Intercept,diagnosis1) 9190
## cor(Intercept,diagnosis2) 12435
## cor(diagnosis1,diagnosis2) 12092
##
## Regression Coefficients:
## Estimate Est.Error l-95% CI u-95% CI
## Intercept 6.17 0.03 6.12 6.22
## diagnosis1 0.02 0.03 -0.03 0.07
## diagnosis2 0.02 0.03 -0.03 0.07
## expected1 -0.04 0.01 -0.05 -0.03
## phase1 0.02 0.01 -0.01 0.04
## phase2 -0.01 0.01 -0.02 0.01
## difficulty1 -0.03 0.01 -0.05 -0.02
## difficulty2 -0.00 0.01 -0.02 0.01
## diagnosis1:expected1 -0.01 0.01 -0.02 0.00
## diagnosis2:expected1 -0.00 0.00 -0.01 0.01
## diagnosis1:phase1 -0.01 0.01 -0.03 0.01
## diagnosis2:phase1 -0.00 0.01 -0.02 0.02
## diagnosis1:phase2 -0.01 0.01 -0.02 0.01
## diagnosis2:phase2 0.00 0.01 -0.01 0.02
## expected1:phase1 -0.04 0.01 -0.06 -0.02
## expected1:phase2 0.00 0.01 -0.01 0.02
## diagnosis1:difficulty1 -0.01 0.01 -0.02 0.00
## diagnosis2:difficulty1 0.00 0.01 -0.01 0.02
## diagnosis1:difficulty2 0.02 0.01 0.01 0.03
## diagnosis2:difficulty2 -0.02 0.01 -0.03 -0.01
## expected1:difficulty1 -0.01 0.01 -0.03 0.00
## expected1:difficulty2 0.00 0.01 -0.01 0.02
## phase1:difficulty1 -0.01 0.01 -0.03 0.02
## phase2:difficulty1 -0.00 0.01 -0.02 0.02
## phase1:difficulty2 0.01 0.01 -0.01 0.04
## phase2:difficulty2 -0.01 0.01 -0.03 0.01
## diagnosis1:expected1:phase1 -0.00 0.01 -0.02 0.01
## diagnosis2:expected1:phase1 0.01 0.01 -0.00 0.02
## diagnosis1:expected1:phase2 0.00 0.01 -0.01 0.01
## diagnosis2:expected1:phase2 -0.00 0.01 -0.01 0.01
## diagnosis1:expected1:difficulty1 -0.01 0.01 -0.02 0.01
## diagnosis2:expected1:difficulty1 0.01 0.01 -0.00 0.02
## diagnosis1:expected1:difficulty2 -0.00 0.01 -0.01 0.01
## diagnosis2:expected1:difficulty2 0.00 0.01 -0.01 0.01
## diagnosis1:phase1:difficulty1 0.00 0.01 -0.01 0.02
## diagnosis2:phase1:difficulty1 0.00 0.01 -0.01 0.02
## diagnosis1:phase2:difficulty1 -0.01 0.01 -0.02 0.00
## diagnosis2:phase2:difficulty1 -0.00 0.01 -0.02 0.01
## diagnosis1:phase1:difficulty2 0.00 0.01 -0.01 0.02
## diagnosis2:phase1:difficulty2 -0.01 0.01 -0.03 0.00
## diagnosis1:phase2:difficulty2 -0.00 0.01 -0.02 0.01
## diagnosis2:phase2:difficulty2 0.01 0.01 -0.00 0.03
## expected1:phase1:difficulty1 -0.01 0.01 -0.03 0.02
## expected1:phase2:difficulty1 -0.00 0.01 -0.02 0.02
## expected1:phase1:difficulty2 0.01 0.01 -0.01 0.03
## expected1:phase2:difficulty2 -0.01 0.01 -0.03 0.01
## diagnosis1:expected1:phase1:difficulty1 0.00 0.01 -0.01 0.02
## diagnosis2:expected1:phase1:difficulty1 0.00 0.01 -0.01 0.02
## diagnosis1:expected1:phase2:difficulty1 0.01 0.01 -0.00 0.02
## diagnosis2:expected1:phase2:difficulty1 -0.01 0.01 -0.02 0.00
## diagnosis1:expected1:phase1:difficulty2 0.00 0.01 -0.01 0.02
## diagnosis2:expected1:phase1:difficulty2 -0.00 0.01 -0.02 0.01
## diagnosis1:expected1:phase2:difficulty2 -0.01 0.01 -0.02 0.00
## diagnosis2:expected1:phase2:difficulty2 0.00 0.01 -0.01 0.01
## Rhat Bulk_ESS Tail_ESS
## Intercept 1.00 1257 2922
## diagnosis1 1.00 1384 3130
## diagnosis2 1.00 1632 3199
## expected1 1.00 5728 8644
## phase1 1.00 5107 8540
## phase2 1.00 5216 9337
## difficulty1 1.00 4814 8336
## difficulty2 1.00 4858 8799
## diagnosis1:expected1 1.00 10348 13271
## diagnosis2:expected1 1.00 10286 12734
## diagnosis1:phase1 1.00 5440 9226
## diagnosis2:phase1 1.00 5033 8789
## diagnosis1:phase2 1.00 9242 11314
## diagnosis2:phase2 1.00 9484 11579
## expected1:phase1 1.00 5887 9381
## expected1:phase2 1.00 5239 9521
## diagnosis1:difficulty1 1.00 9997 11062
## diagnosis2:difficulty1 1.00 10562 12419
## diagnosis1:difficulty2 1.00 11211 12983
## diagnosis2:difficulty2 1.00 11525 13218
## expected1:difficulty1 1.00 4894 8458
## expected1:difficulty2 1.00 5300 8739
## phase1:difficulty1 1.00 4742 8273
## phase2:difficulty1 1.00 4199 8491
## phase1:difficulty2 1.00 4816 9280
## phase2:difficulty2 1.00 4547 7937
## diagnosis1:expected1:phase1 1.00 13718 14036
## diagnosis2:expected1:phase1 1.00 13722 13021
## diagnosis1:expected1:phase2 1.00 10816 12500
## diagnosis2:expected1:phase2 1.00 11218 13346
## diagnosis1:expected1:difficulty1 1.00 10306 13500
## diagnosis2:expected1:difficulty1 1.00 11979 13393
## diagnosis1:expected1:difficulty2 1.00 11054 12978
## diagnosis2:expected1:difficulty2 1.00 11268 12715
## diagnosis1:phase1:difficulty1 1.00 9380 12504
## diagnosis2:phase1:difficulty1 1.00 10149 12743
## diagnosis1:phase2:difficulty1 1.00 9314 11792
## diagnosis2:phase2:difficulty1 1.00 9906 12839
## diagnosis1:phase1:difficulty2 1.00 9038 12396
## diagnosis2:phase1:difficulty2 1.00 9988 13029
## diagnosis1:phase2:difficulty2 1.00 9117 13078
## diagnosis2:phase2:difficulty2 1.00 9025 12760
## expected1:phase1:difficulty1 1.00 5180 8462
## expected1:phase2:difficulty1 1.00 4292 7684
## expected1:phase1:difficulty2 1.00 5284 8058
## expected1:phase2:difficulty2 1.00 4688 7872
## diagnosis1:expected1:phase1:difficulty1 1.00 9765 13094
## diagnosis2:expected1:phase1:difficulty1 1.00 10073 12791
## diagnosis1:expected1:phase2:difficulty1 1.00 9494 12343
## diagnosis2:expected1:phase2:difficulty1 1.00 9763 12819
## diagnosis1:expected1:phase1:difficulty2 1.00 9227 12066
## diagnosis2:expected1:phase1:difficulty2 1.00 9940 12701
## diagnosis1:expected1:phase2:difficulty2 1.00 9034 13206
## diagnosis2:expected1:phase2:difficulty2 1.00 9861 13223
##
## Further Distributional Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 0.25 0.00 0.24 0.25 1.00 11159 11907
## ndt 106.89 5.16 96.49 116.71 1.00 10792 12313
##
## Draws were sampled using sample(hmc). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).# get the estimates and compute group comparisons
df.m.pal = post.draws %>%
select(starts_with("b_")) %>%
mutate(
b_postvolatile = - b_phase1 - b_phase2,
b_difficult = - b_difficulty1 - b_difficulty2
)
# plot the posterior distributions
df.m.pal %>%
pivot_longer(cols = starts_with("b_"), names_to = "coef", values_to = "estimate") %>%
subset(!startsWith(coef, "b_Int")) %>%
mutate(
coef = substr(coef, 3, nchar(coef)),
coef = str_replace_all(coef, ":", " x "),
coef = str_replace_all(coef, "diagnosis1", "ADHD"),
coef = str_replace_all(coef, "diagnosis2", "ADHD+ASD"),
coef = str_replace_all(coef, "expected1", "expected"),
coef = str_replace_all(coef, "expected2", "unexpected"),
coef = str_replace_all(coef, "phase1", "prevolatile"),
coef = str_replace_all(coef, "phase2", "volatile"),
coef = str_replace_all(coef, "difficulty1", "easy"),
coef = str_replace_all(coef, "difficulty2", "medium"),
coef = fct_reorder(coef, desc(estimate))
) %>%
group_by(coef) %>%
mutate(
cred = case_when(
(mean(estimate) < 0 & quantile(estimate, probs = 0.975) < 0) |
(mean(estimate) > 0 & quantile(estimate, probs = 0.025) > 0) ~ "credible",
T ~ "not credible"
)
) %>% ungroup() %>%
ggplot(aes(x = estimate, y = coef, fill = cred)) +
geom_vline(xintercept = 0, linetype = 'dashed') +
ggdist::stat_halfeye(alpha = 0.7) + ylab(NULL) + theme_bw() +
scale_fill_manual(values = c(c_dark, c_light)) + theme(legend.position = "none")
## create design matrix to figure out how to set comparisons for hypotheses
df.des = cbind(df.pal %>% select(diagnosis, expected, phase, difficulty),
model.matrix(~ diagnosis * expected * phase * difficulty,
data = df.pal)) %>% distinct()
# H1c: COMP(unexp-exp) != ADHD(unexp-exp)
as.data.frame(t(df.des %>%
filter(diagnosis != "BOTH") %>%
group_by(diagnosis, expected) %>%
summarise(across(where(is.numeric), ~ mean(.x))) %>%
group_by(diagnosis) %>%
summarise(across(where(is.numeric), ~ diff(.x))) %>% # unexpected - expected
select(where(is.numeric)) %>%
map_df(~ diff(.x)))) %>%
filter(V1 != 0) # COMP - ADHD## V1
## diagnosis1:expected1 4
## diagnosis2:expected1 2h1c = hypothesis(m.pal, "0 > 2*diagnosis1:expected1 + diagnosis2:expected1")
h1c## Hypothesis Tests for class b:
## Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio
## 1 (0)-(2*diagnosis1... > 0 0.01 0.01 0 0.03 9.44
## Post.Prob Star
## 1 0.9
## ---
## 'CI': 90%-CI for one-sided and 95%-CI for two-sided hypotheses.
## '*': For one-sided hypotheses, the posterior probability exceeds 95%;
## for two-sided hypotheses, the value tested against lies outside the 95%-CI.
## Posterior probabilities of point hypotheses assume equal prior probabilities.# COMP(unexp-exp) != BOTH(unexp-exp)
e.both = hypothesis(m.pal, "0 > diagnosis1:expected1 + 2*diagnosis2:expected1", alpha = 0.025)
e.both## Hypothesis Tests for class b:
## Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio
## 1 (0)-(diagnosis1:e... > 0 0.01 0.01 -0.01 0.02 3.49
## Post.Prob Star
## 1 0.78
## ---
## 'CI': 95%-CI for one-sided and 97.5%-CI for two-sided hypotheses.
## '*': For one-sided hypotheses, the posterior probability exceeds 97.5%;
## for two-sided hypotheses, the value tested against lies outside the 97.5%-CI.
## Posterior probabilities of point hypotheses assume equal prior probabilities.# ADHD(unexp-exp) != BOTH(unexp-exp)
e.clin = hypothesis(m.pal, "diagnosis1:expected1 < diagnosis2:expected1", alpha = 0.025)
e.clin## Hypothesis Tests for class b:
## Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio
## 1 (diagnosis1:expec... < 0 0 0.01 -0.02 0.01 2.5
## Post.Prob Star
## 1 0.71
## ---
## 'CI': 95%-CI for one-sided and 97.5%-CI for two-sided hypotheses.
## '*': For one-sided hypotheses, the posterior probability exceeds 97.5%;
## for two-sided hypotheses, the value tested against lies outside the 97.5%-CI.
## Posterior probabilities of point hypotheses assume equal prior probabilities.## exploration: task effects
# expected versus unexpected
e1 = hypothesis(m.pal, "0 > 2*expected1", alpha = 0.025)
e1## Hypothesis Tests for class b:
## Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio
## 1 (0)-(2*expected1) > 0 0.08 0.01 0.06 0.11 Inf
## Post.Prob Star
## 1 1 *
## ---
## 'CI': 95%-CI for one-sided and 97.5%-CI for two-sided hypotheses.
## '*': For one-sided hypotheses, the posterior probability exceeds 97.5%;
## for two-sided hypotheses, the value tested against lies outside the 97.5%-CI.
## Posterior probabilities of point hypotheses assume equal prior probabilities.# volatile versus prevolatile
e2 = hypothesis(m.pal, "0 < phase1 - phase2", alpha = 0.025)
e2## Hypothesis Tests for class b:
## Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio
## 1 (0)-(phase1-phase2) < 0 -0.02 0.02 -0.06 0.01 15.03
## Post.Prob Star
## 1 0.94
## ---
## 'CI': 95%-CI for one-sided and 97.5%-CI for two-sided hypotheses.
## '*': For one-sided hypotheses, the posterior probability exceeds 97.5%;
## for two-sided hypotheses, the value tested against lies outside the 97.5%-CI.
## Posterior probabilities of point hypotheses assume equal prior probabilities.# volatile versus postvolatile
e3 = hypothesis(m.pal, "0 < phase1 + 2*phase2", alpha = 0.025)
e3## Hypothesis Tests for class b:
## Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio
## 1 (0)-(phase1+2*pha... < 0 0 0.02 -0.03 0.03 1.15
## Post.Prob Star
## 1 0.54
## ---
## 'CI': 95%-CI for one-sided and 97.5%-CI for two-sided hypotheses.
## '*': For one-sided hypotheses, the posterior probability exceeds 97.5%;
## for two-sided hypotheses, the value tested against lies outside the 97.5%-CI.
## Posterior probabilities of point hypotheses assume equal prior probabilities.# difficult versus medium
e4 = hypothesis(m.pal, "0 < -difficulty1 - 2*difficulty2", alpha = 0.025)
e4## Hypothesis Tests for class b:
## Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio
## 1 (0)-(-difficulty1... < 0 -0.04 0.02 -0.07 -0.01 366.35
## Post.Prob Star
## 1 1 *
## ---
## 'CI': 95%-CI for one-sided and 97.5%-CI for two-sided hypotheses.
## '*': For one-sided hypotheses, the posterior probability exceeds 97.5%;
## for two-sided hypotheses, the value tested against lies outside the 97.5%-CI.
## Posterior probabilities of point hypotheses assume equal prior probabilities.# medium versus easy
e5 = hypothesis(m.pal, "0 < -difficulty1 + difficulty2", alpha = 0.025)
e5## Hypothesis Tests for class b:
## Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio
## 1 (0)-(-difficulty1... < 0 -0.03 0.01 -0.06 0 44.92
## Post.Prob Star
## 1 0.98 *
## ---
## 'CI': 95%-CI for one-sided and 97.5%-CI for two-sided hypotheses.
## '*': For one-sided hypotheses, the posterior probability exceeds 97.5%;
## for two-sided hypotheses, the value tested against lies outside the 97.5%-CI.
## Posterior probabilities of point hypotheses assume equal prior probabilities.# difficult versus easy
e6 = hypothesis(m.pal, "0 < -2*difficulty1 - difficulty2", alpha = 0.025)
e6## Hypothesis Tests for class b:
## Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio
## 1 (0)-(-2*difficult... < 0 -0.07 0.02 -0.1 -0.04 Inf
## Post.Prob Star
## 1 1 *
## ---
## 'CI': 95%-CI for one-sided and 97.5%-CI for two-sided hypotheses.
## '*': For one-sided hypotheses, the posterior probability exceeds 97.5%;
## for two-sided hypotheses, the value tested against lies outside the 97.5%-CI.
## Posterior probabilities of point hypotheses assume equal prior probabilities.## extract predicted differences in ms instead of log data
df.new = df.pal %>%
select(diagnosis, phase, expected, difficulty) %>%
distinct() %>%
mutate(
condition = paste(diagnosis, phase, expected, difficulty, sep = "_")
)
df.ms = as.data.frame(
fitted(m.pal, summary = F,
newdata = df.new %>% select(diagnosis, phase, expected, difficulty),
re_formula = NA))
colnames(df.ms) = df.new$condition
# calculate our difference columns
df.ms = df.ms %>%
mutate(
COMP_expected = rowMeans(across(matches("COMP_.*_expected_.*"))),
COMP_unexpected = rowMeans(across(matches("COMP_.*_unexpected_.*"))),
ADHD_expected = rowMeans(across(matches("ADHD_.*_expected_.*"))),
ADHD_unexpected = rowMeans(across(matches("ADHD_.*_unexpected_.*"))),
BOTH_expected = rowMeans(across(matches("BOTH_.*_expected_.*"))),
BOTH_unexpected = rowMeans(across(matches("BOTH_.*_unexpected_.*"))),
COMP_unexp_exp = COMP_unexpected - COMP_expected,
ADHD_unexp_exp = ADHD_unexpected - ADHD_expected,
BOTH_unexp_exp = BOTH_unexpected - BOTH_expected,
`h1c_COMP-ADHD unexpected-expected` = ADHD_unexp_exp - COMP_unexp_exp,
`e_COMP-BOTH unexpected-expected` = BOTH_unexp_exp - COMP_unexp_exp,
`e_ADHD-BOTH unexpected-expected` = ADHD_unexp_exp - BOTH_unexp_exp,
`e1_unexpected-expected` = rowMeans(across(matches(".*_unexpected_.*"))) -
rowMeans(across(matches(".*_expected_.*"))),
`e2_volatile-prevolatile` = rowMeans(across(matches(".*_volatile_.*"))) -
rowMeans(across(matches(".*_prevolatile_.*"))),
`e3_volatile-postvolatile` = rowMeans(across(matches(".*_volatile_.*"))) -
rowMeans(across(matches(".*_postvolatile_.*"))),
`e4_difficult-medium` = rowMeans(across(matches(".*_difficult"))) -
rowMeans(across(matches(".*_medium"))),
`e5_medium-easy` = rowMeans(across(matches(".*_medium"))) -
rowMeans(across(matches(".*_easy"))),
`e6_difficult-easy` = rowMeans(across(matches(".*_difficult"))) -
rowMeans(across(matches(".*_easy")))
)
kable(df.ms %>%
select(starts_with("e") | starts_with("h")) %>%
summarise_all(list(lower.ci = lower_ci, mean = mean, upper.ci = upper_ci)) %>%
t %>%
as.data.frame %>%
rownames_to_column() %>%
separate(rowname, into = c("id", "comparison", "stat"), sep = "_") %>%
pivot_wider(names_from = stat, values_from = V1) %>%
arrange(comparison), digits = 2)| id | comparison | lower.ci | mean | upper.ci |
|---|---|---|---|---|
| e | ADHD-BOTH unexpected-expected | -12.73 | 4.50 | 21.66 |
| h1c | COMP-ADHD unexpected-expected | -4.37 | 12.89 | 30.22 |
| e | COMP-BOTH unexpected-expected | -8.61 | 8.39 | 25.57 |
| e6 | difficult-easy | -11.55 | 3.47 | 18.97 |
| e4 | difficult-medium | 5.59 | 20.76 | 36.26 |
| e5 | medium-easy | -0.26 | 14.23 | 28.62 |
| e1 | unexpected-expected | 27.75 | 40.71 | 54.22 |
| e3 | volatile-postvolatile | -14.98 | 0.80 | 16.60 |
| e2 | volatile-prevolatile | -29.65 | -13.41 | 2.76 |
equivalence_test(df.ms %>% select(`h1c_COMP-ADHD unexpected-expected`,
`e_COMP-BOTH unexpected-expected`,
`e_ADHD-BOTH unexpected-expected`),
range = rope_range(m.pal))## # Test for Practical Equivalence
##
## ROPE: [-17.71 17.71]
##
## Parameter | H0 | inside ROPE | 95% HDI
## -----------------------------------------------------------------------------
## h1c_COMP-ADHD unexpected-expected | Undecided | 71.64 % | [-4.37, 30.22]
## e_COMP-BOTH unexpected-expected | Undecided | 88.00 % | [-8.61, 25.57]
## e_ADHD-BOTH unexpected-expected | Undecided | 95.87 % | [-12.73, 21.66]# calculate effect sizes
df.effect = post.draws %>%
mutate(across(starts_with("sd")|starts_with("sigma"), ~.^2)) %>%
mutate(
sumvar = sqrt(rowSums(select(., starts_with("sd")|starts_with("sigma")))),
h1c = -(2*`b_diagnosis1:expected1` + `b_diagnosis2:expected1`) / sumvar,
e1 = -(2*b_expected1) / sumvar,
e6 = -(2*`b_difficulty1` + `b_difficulty2`) / sumvar
)
kable(df.effect %>% select(starts_with("e")|starts_with("h")) %>%
pivot_longer(cols = everything(), values_to = "estimate") %>%
group_by(name) %>%
summarise(
ci.lo = lower_ci(estimate),
mean = mean(estimate),
ci.hi = upper_ci(estimate),
interpret = interpret_cohens_d(mean)
), digits = 3
)| name | ci.lo | mean | ci.hi | interpret |
|---|---|---|---|---|
| e1 | 0.169 | 0.249 | 0.331 | small |
| e6 | 0.130 | 0.219 | 0.312 | small |
| h1c | -0.018 | 0.035 | 0.089 | very small |
h1c: estimate = 0.01 [0, 0.03], posterior probability = 90.42%
3.3 Plots
# rain cloud plot including all factors
df.pal %>%
ggplot(aes(diagnosis, rt.cor, fill = expected, colour = expected)) + #
geom_rain(rain.side = 'r',
boxplot.args = list(color = "black", outlier.shape = NA, show.legend = FALSE, alpha = .8),
violin.args = list(color = "black", outlier.shape = NA, alpha = .8),
boxplot.args.pos = list(
position = ggpp::position_dodgenudge(x = 0, width = 0.3), width = 0.3
),
point.args = list(show_guide = FALSE, alpha = .5),
violin.args.pos = list(
width = 0.6, position = position_nudge(x = 0.16)),
point.args.pos = list(position = ggpp::position_dodgenudge(x = -0.25, width = 0.1))) +
scale_fill_manual(values = col.cnd) +
scale_color_manual(values = col.cnd) +
#ylim(0, 1) +
labs(title = "Reaction times per subject", x = "", y = "rt (ms)") +
facet_wrap(difficulty ~ phase) +
theme_bw() +
theme(legend.position = "bottom", plot.title = element_text(hjust = 0.5),
legend.direction = "horizontal", text = element_text(size = 15))
# rain cloud plot excluding difficulty
df.pal %>%
group_by(subID, diagnosis, expected, phase) %>%
summarise(
rt.cor = median(rt.cor, na.rm = T)
) %>%
group_by(subID, diagnosis, phase) %>%
summarise(
rt.exp = diff(rt.cor)
) %>%
ggplot(aes(phase, rt.exp, fill = diagnosis, colour = diagnosis)) + #
geom_rain(rain.side = 'r',
boxplot.args = list(color = "black", outlier.shape = NA, show.legend = FALSE, alpha = .8),
violin.args = list(color = "black", outlier.shape = NA, alpha = .8),
boxplot.args.pos = list(
position = ggpp::position_dodgenudge(x = 0, width = 0.3), width = 0.3
),
point.args = list(show_guide = FALSE, alpha = .5),
violin.args.pos = list(
width = 0.6, position = position_nudge(x = 0.16)),
point.args.pos = list(position = ggpp::position_dodgenudge(x = -0.25, width = 0.1))) +
scale_fill_manual(values = col.grp) +
scale_color_manual(values = col.grp) +
#ylim(0, 1) +
labs(title = "", x = "", y = "rt difference (ms)") +
theme_bw() +
theme(legend.position = "bottom", plot.title = element_blank(),
legend.direction = "horizontal", text = element_text(size = 15),
legend.title = element_blank())
ggsave("plots/FigRT.svg", units = "cm", width = 27, height = 9)
# plot transition effects specifically
df.pal %>%
group_by(subID, diagnosis, phase, expected) %>%
summarise(rt.cor = median(rt.cor, na.rm = T)) %>%
group_by(diagnosis, phase, expected) %>%
summarise(
rt.mn = mean(rt.cor),
rt.se = sd(rt.cor)
) %>%
ggplot(aes(y = rt.mn, x = phase,
group = diagnosis, colour = diagnosis)) +
geom_line(position = position_dodge(0.4), linewidth = 1) +
geom_errorbar(aes(ymin = rt.mn - rt.se,
ymax = rt.mn + rt.se), linewidth = 1, width = 0.5,
position = position_dodge(0.4)) +
geom_point(position = position_dodge(0.4), size = 5) +
labs (x = "phase", y = "rt (ms)") +
ggtitle("Diagnosis x phase x expected") +
facet_wrap(. ~ expected) +
scale_fill_manual(values = col.grp) +
scale_color_manual(values = col.grp) +
theme_bw() +
theme(legend.position = "bottom", plot.title = element_text(hjust = 0.5))
4 Explorative analysis of errors
Last but not least, we are going to explore possible differences with regards to mean accuracies using a Bayesian ANOVA.
# create accuracies dataframe
df.acc = df.tsk %>% filter(expected != "neutral") %>% droplevels() %>%
group_by(subID, diagnosis, phase, expected, difficulty) %>%
summarise(acc = mean(acc)*100)
# check normal distribution
kable(df.tsk %>% group_by(subID, diagnosis, phase, expected, difficulty) %>%
summarise(acc = mean(acc)*100) %>%
group_by(diagnosis, phase, expected, difficulty) %>%
shapiro_test(acc) %>%
mutate(
sig = if_else(p < 0.05, "*", "")
))| diagnosis | phase | expected | difficulty | variable | statistic | p | sig |
|---|---|---|---|---|---|---|---|
| ADHD | prevolatile | expected | easy | acc | 0.7530551 | 0.0001005 | * |
| ADHD | prevolatile | expected | medium | acc | 0.7069895 | 0.0000240 | * |
| ADHD | prevolatile | expected | difficult | acc | 0.8478059 | 0.0031053 | * |
| ADHD | prevolatile | unexpected | easy | acc | 0.7684096 | 0.0001668 | * |
| ADHD | prevolatile | unexpected | medium | acc | 0.6030083 | 0.0000014 | * |
| ADHD | prevolatile | unexpected | difficult | acc | 0.6419123 | 0.0000039 | * |
| ADHD | volatile | expected | easy | acc | 0.7936914 | 0.0003993 | * |
| ADHD | volatile | expected | medium | acc | 0.8355305 | 0.0019035 | * |
| ADHD | volatile | expected | difficult | acc | 0.8366283 | 0.0019874 | * |
| ADHD | volatile | neutral | easy | acc | 0.7837217 | 0.0002813 | * |
| ADHD | volatile | neutral | medium | acc | 0.8443747 | 0.0027041 | * |
| ADHD | volatile | neutral | difficult | acc | 0.8929161 | 0.0215045 | * |
| ADHD | volatile | unexpected | easy | acc | 0.7557361 | 0.0001096 | * |
| ADHD | volatile | unexpected | medium | acc | 0.8294574 | 0.0015023 | * |
| ADHD | volatile | unexpected | difficult | acc | 0.8536698 | 0.0039446 | * |
| ADHD | postvolatile | expected | easy | acc | 0.8440215 | 0.0026661 | * |
| ADHD | postvolatile | expected | medium | acc | 0.9104130 | 0.0482689 | * |
| ADHD | postvolatile | expected | difficult | acc | 0.8661620 | 0.0066462 | * |
| ADHD | postvolatile | unexpected | easy | acc | 0.5959201 | 0.0000012 | * |
| ADHD | postvolatile | unexpected | medium | acc | 0.5555850 | 0.0000005 | * |
| ADHD | postvolatile | unexpected | difficult | acc | 0.7380254 | 0.0000621 | * |
| BOTH | prevolatile | expected | easy | acc | 0.7693266 | 0.0001720 | * |
| BOTH | prevolatile | expected | medium | acc | 0.8136402 | 0.0008245 | * |
| BOTH | prevolatile | expected | difficult | acc | 0.8606902 | 0.0052777 | * |
| BOTH | prevolatile | unexpected | easy | acc | 0.6126826 | 0.0000018 | * |
| BOTH | prevolatile | unexpected | medium | acc | 0.5222700 | 0.0000002 | * |
| BOTH | prevolatile | unexpected | difficult | acc | 0.8256146 | 0.0012958 | * |
| BOTH | volatile | expected | easy | acc | 0.8657282 | 0.0065251 | * |
| BOTH | volatile | expected | medium | acc | 0.7762210 | 0.0002173 | * |
| BOTH | volatile | expected | difficult | acc | 0.9051781 | 0.0377818 | * |
| BOTH | volatile | neutral | easy | acc | 0.7245658 | 0.0000409 | * |
| BOTH | volatile | neutral | medium | acc | 0.6932121 | 0.0000161 | * |
| BOTH | volatile | neutral | difficult | acc | 0.9250169 | 0.0966207 | |
| BOTH | volatile | unexpected | easy | acc | 0.6507543 | 0.0000049 | * |
| BOTH | volatile | unexpected | medium | acc | 0.7485284 | 0.0000868 | * |
| BOTH | volatile | unexpected | difficult | acc | 0.7445783 | 0.0000765 | * |
| BOTH | postvolatile | expected | easy | acc | 0.7509063 | 0.0000937 | * |
| BOTH | postvolatile | expected | medium | acc | 0.8264958 | 0.0013403 | * |
| BOTH | postvolatile | expected | difficult | acc | 0.9085878 | 0.0443054 | * |
| BOTH | postvolatile | unexpected | easy | acc | 0.4955981 | 0.0000001 | * |
| BOTH | postvolatile | unexpected | medium | acc | 0.5505637 | 0.0000004 | * |
| BOTH | postvolatile | unexpected | difficult | acc | 0.7628110 | 0.0001384 | * |
| COMP | prevolatile | expected | easy | acc | 0.6048243 | 0.0000015 | * |
| COMP | prevolatile | expected | medium | acc | 0.6717248 | 0.0000087 | * |
| COMP | prevolatile | expected | difficult | acc | 0.9220654 | 0.0838957 | |
| COMP | prevolatile | unexpected | easy | acc | 0.5904384 | 0.0000010 | * |
| COMP | prevolatile | unexpected | medium | acc | 0.5555850 | 0.0000005 | * |
| COMP | prevolatile | unexpected | difficult | acc | 0.8469958 | 0.0030052 | * |
| COMP | volatile | expected | easy | acc | 0.8053882 | 0.0006083 | * |
| COMP | volatile | expected | medium | acc | 0.7784453 | 0.0002345 | * |
| COMP | volatile | expected | difficult | acc | 0.8934097 | 0.0219914 | * |
| COMP | volatile | neutral | easy | acc | 0.6048243 | 0.0000015 | * |
| COMP | volatile | neutral | medium | acc | 0.7416438 | 0.0000696 | * |
| COMP | volatile | neutral | difficult | acc | 0.8721173 | 0.0085735 | * |
| COMP | volatile | unexpected | easy | acc | 0.6046337 | 0.0000015 | * |
| COMP | volatile | unexpected | medium | acc | 0.6449235 | 0.0000042 | * |
| COMP | volatile | unexpected | difficult | acc | 0.7733926 | 0.0001974 | * |
| COMP | postvolatile | expected | easy | acc | 0.6488874 | 0.0000047 | * |
| COMP | postvolatile | expected | medium | acc | 0.7015343 | 0.0000205 | * |
| COMP | postvolatile | expected | difficult | acc | 0.8164406 | 0.0009154 | * |
| COMP | postvolatile | unexpected | easy | acc | 0.4735522 | 0.0000001 | * |
| COMP | postvolatile | unexpected | medium | acc | 0.4225841 | 0.0000000 | * |
| COMP | postvolatile | unexpected | difficult | acc | 0.7202060 | 0.0000358 | * |
# rank transform the data
df.acc = df.acc %>%
ungroup() %>%
mutate(
racc = rank(acc)
)
# run the ANOVA
if (!file.exists(file.path(brms_dir, "aov_acc.rds"))){
aov.acc = anovaBF(racc ~ diagnosis * phase * expected * difficulty, data = df.acc)
saveRDS(aov.acc, file = file.path(brms_dir, "aov_acc.rds"))
} else {
aov.acc = readRDS(file.path(brms_dir, "aov_acc.rds"))
}
bf.acc = extractBF(aov.acc, logbf = T)
kable(head(bf.acc %>% arrange(desc(bf))))| bf | error | time | code | |
|---|---|---|---|---|
| diagnosis + expected + difficulty | 46.41268 | 0.0176216 | Fri Sep 5 15:14:41 2025 | 23cdf312fa7 |
| diagnosis + difficulty | 46.24773 | 0.0220083 | Fri Sep 5 15:14:41 2025 | 23cd20502937 |
| diagnosis + expected + difficulty + diagnosis:difficulty | 45.01890 | 0.0189216 | Fri Sep 5 15:14:44 2025 | 23cd306bedbc |
| diagnosis + difficulty + diagnosis:difficulty | 44.85054 | 0.0153529 | Fri Sep 5 15:14:43 2025 | 23cd186e77cb |
| diagnosis + phase + expected + phase:expected + difficulty | 44.76534 | 0.1771458 | Fri Sep 5 15:14:42 2025 | 23cd2da83298 |
| expected + difficulty | 44.41762 | 0.0165907 | Fri Sep 5 15:14:41 2025 | 23cd4f1e8132 |
# print overall accuracy rates for all the effects included in the best model
df.acc %>%
group_by(diagnosis, difficulty, expected) %>%
summarise(mean_accuracy = mean(acc, na.rm = T),
sd_accuracy = sd(acc, na.rm = T))## # A tibble: 18 × 5
## # Groups: diagnosis, difficulty [9]
## diagnosis difficulty expected mean_accuracy sd_accuracy
## <fct> <fct> <fct> <dbl> <dbl>
## 1 ADHD easy expected 94.1 6.87
## 2 ADHD easy unexpected 88.4 17.0
## 3 ADHD medium expected 92.6 8.14
## 4 ADHD medium unexpected 90.7 14.3
## 5 ADHD difficult expected 88.6 9.62
## 6 ADHD difficult unexpected 86.6 16.0
## 7 BOTH easy expected 94.9 6.58
## 8 BOTH easy unexpected 92.4 12.8
## 9 BOTH medium expected 93.9 7.07
## 10 BOTH medium unexpected 91.9 13.8
## 11 BOTH difficult expected 87.9 9.95
## 12 BOTH difficult unexpected 83.0 18.6
## 13 COMP easy expected 97.9 3.23
## 14 COMP easy unexpected 93.9 10.8
## 15 COMP medium expected 96.3 5.34
## 16 COMP medium unexpected 94.7 11.4
## 17 COMP difficult expected 89.8 8.68
## 18 COMP difficult unexpected 84.8 16.6df.acc %>%
group_by(diagnosis) %>%
summarise(mean_accuracy = mean(acc, na.rm = T),
sd_accuracy = sd(acc, na.rm = T))## # A tibble: 3 × 3
## diagnosis mean_accuracy sd_accuracy
## <fct> <dbl> <dbl>
## 1 ADHD 90.2 12.8
## 2 BOTH 90.7 12.8
## 3 COMP 92.9 11.1df.acc %>%
group_by(expected) %>%
summarise(mean_accuracy = mean(acc, na.rm = T),
sd_accuracy = sd(acc, na.rm = T))## # A tibble: 2 × 3
## expected mean_accuracy sd_accuracy
## <fct> <dbl> <dbl>
## 1 expected 92.9 8.18
## 2 unexpected 89.6 15.2df.acc %>%
group_by(difficulty) %>%
summarise(mean_accuracy = mean(acc, na.rm = T),
sd_accuracy = sd(acc, na.rm = T))## # A tibble: 3 × 3
## difficulty mean_accuracy sd_accuracy
## <fct> <dbl> <dbl>
## 1 easy 93.6 10.9
## 2 medium 93.4 10.7
## 3 difficult 86.8 13.9Accuracies were generally high, with a grand average of 91.25% accurate responses across diagnostic groups. Accuracies seems to have differed between diagnostic groups. Let’s do some plotting to figure out where our differences lie.
4.1 Plots
# rain cloud plot including all factors
df.acc %>%
ggplot(aes(expected, acc, fill = diagnosis, colour = diagnosis)) + #
geom_rain(rain.side = 'r',
boxplot.args = list(color = "black", outlier.shape = NA, show.legend = FALSE, alpha = .8),
violin.args = list(color = "black", outlier.shape = NA, alpha = .8),
boxplot.args.pos = list(
position = ggpp::position_dodgenudge(x = 0, width = 0.3), width = 0.3
),
point.args = list(show_guide = FALSE, alpha = .5),
violin.args.pos = list(
width = 0.6, position = position_nudge(x = 0.16)),
point.args.pos = list(position = ggpp::position_dodgenudge(x = -0.25, width = 0.1))) +
scale_fill_manual(values = col.grp) +
scale_color_manual(values = col.grp) +
labs(title = "Accuracies per subject", x = "", y = "accuracy (%)") +
facet_wrap(phase ~ difficulty) +
theme_bw() +
theme(legend.position = "bottom", plot.title = element_text(hjust = 0.5),
legend.direction = "horizontal", text = element_text(size = 15))
# plot as line plot without the phases
df.acc %>%
group_by(subID, diagnosis, expected, difficulty) %>%
summarise(
acc = mean(acc, na.rm = T)
) %>%
group_by(diagnosis, expected, difficulty) %>%
summarise(
acc.mn = mean(acc),
acc.se = sd(acc)
) %>%
ggplot(aes(y = acc.mn, x = difficulty, group = expected, colour = expected)) +
geom_line(position = position_dodge(0.4), linewidth = 1) +
geom_errorbar(aes(ymin = acc.mn - acc.se,
ymax = acc.mn + acc.se), linewidth = 1, width = 0.5,
position = position_dodge(0.4)) +
geom_point(position = position_dodge(0.4), size = 5) +
facet_wrap(. ~ diagnosis) +
#ylim(75,100) +
labs (x = "difficulty", y = "accuracy") +
scale_fill_manual(values = col.cnd) +
scale_color_manual(values = col.cnd) +
theme_bw() +
theme(legend.position = "bottom", plot.title = element_text(hjust = 0.5))